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ABSTRACT 


The ability of a submarine to maintain ordered depth, especially during periscope 
depth operations at low speeds, is vital for the vessel to perform its mission and avoid 
detection. Modern submarines exhibit an inherent phenomenon that produces an 
undesirable ship response at low speeds, commonly referred to as dive plane reversal. The 
physical parameters that govern this occurrence are related in this thesis to the problem of 
multiple steady state solutions in the vertical plane. 

Generic solution branching, in the form of pitchfork bifurcations, can occur when the 
nominal level flight path loses its stability. A systematic study reveals the existence of a 
critical Froude number, based on the vessel's speed and metacentric height, where this 
branching occurs. Bifurcation theory techniques and numerical computations are utilized 
to classify the effect that geometric parameters, trim and ballast conditions, and 


hydrodynamic properties have on the existence of these multiple solutions. 
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I. INTRODUCTION 


One of the most critical functions of a submersible is accurate depth keeping at the 
commanded depth. Such a function can be carried out either manually or automatically, 
especially in cases where human intervention is not possible or desirable. Due to the 
technological significance and numerous scientific applications of submersible vehicle 
systems, design issues of appropriate depth keeping control laws have received wide 
attention in the past. Such control system designs include linear and nonlinear controllers 
[Refs. 1&2], model based compensators [Ref. 3], adaptive control [Ref. 4], and sliding 
mode control laws [Refs. 5&6]. 

Response accuracy and stability are primary considerations in. designing depth 
keeping control law. Of paramount importance, in this area, are the robustness properties 
of the particular design; 1.e., its ability to maintain accuracy and stability in the presence of 
incomplete sensor and environmental oran as well as actual/mathematical model 
mismatch. The scope of the work in this thesis is to demonstrate a potential loss of 
stability that may occur when a submarine is operating at low speeds. This loss of stability 
can occur regardless of the particular means used for depth control. The study is 
accomplished through the use of an eigenvalue analysis, a steady state analysis and a 
controllability analysis [Ref. 7]. It is shown that such a loss of stability is accompanied by a 
slow divergence of trajectories away from the commanded path. Solution branching 
occurs in the form of generic pitchfork bifurcations [Refs. 8, 9]. A complete 
characterization of the problem is given utilizing singularity techniques, which have been 


proven to be very useful in the analysis of similar problems (Refs. 10, 11, 12]. The use of 


bifurcation theory allows the crucial vehicle parameters that govern the problem of 
solution branching to be determined, and the develop guidelines to prevent its occurrence. 
Finallv, a new look at the problem of dive plane reversal [Ref. 13], based on solution 
branching results is presented. The term dive plane reversal refers to a well-known 
phenomenon in submarine operations where, during low speed depth keeping, there is a 
need to reverse the direction of dive plane deflection in order to execute a given change in 
depth. Physically, this can be explained by considering the relative magnitude of the 
hydrodynamic forces. At moderate and high speeds, the normal force on the submarine's 
hull due to the angle of attack exceeds the normal dive plane force and the boat responds 
to ordered dive plane angles as expected. The phenomenon of dive plane reversal occurs 
at speeds below a certain critical speed in. which the normal hull force is less than the 
normal dive plane force and the response of the boat is reversed. [Ref. 14] Vehicle 
modeling in this work follows standard notation [Ref. 15] and numerical results are 
presented for the DARPA SUBOFF model [Ref. 16] for which a set of hydrodynamic 
coefficients and geometric properties is available. Special emphasis is given to identifying 
the proper non-dimensional parameters in the problem, so that extension of these results 
to full scale models and other designs is possible using minimal experimental and/or 


analytical results. 


Il. VEHICLE MODELING 


A. EQUATIONS OF MOTION 
l. Introduction 
For the purpose of this problem, and the subject of the maneuvering and motion 
control of the vehicle, the following assumptions are made: 
1) The vehicle behaves as a rigid body; 
2) The earth's rotation is negligible as far as acceleration components of the 
center of mass are concerned; 
3) The primary forces that act on the vehicle have inertial, gravitational, 
hydrostatic and hydrodynamic origins. 
2. Coordinate Systems And Positional Definitions 
A global navigation frame, OXYZ, às shown in Figure 2.1, is defined with origin, 
O, and a set of axes aligned with directions North, East and Down. This produces a right- 
hand reference frame with unit vectors 7, J, and K. Ignoring the earth's rotation rate in 
comparison to the angular rates produced by the vessel's motion, it can be said that the Да 


J, and K coordinate frame is an inertial reference frame in which Newton's Laws of 


Motion will be valid. As seen in Figure 2.1, a vehicle's position R ‚ їп this frame will have 
the vector components, R.=/X,/ +Y,J+Z,K]. A standard convention that is used 
in both aircraft and marine vehicle dynamics places the Y axis to the right while looking 


along the X axis, and the Z axis is positive downwards. Figure 2.1 also shows a vehicle 


with some general attitude in the original coordinate system. 





Figure 2-1. Coordinate Axes Convention. 


Now define a body fixed frame of reference O 4, with the origin O , located on 
the vehicle centerline, moving and rotating with the vehicle, in which the vehicle's center 
of mass, G, has some position other than the origin of the vehicle fixed frame (refer to 
Figure 2-1). This origin will be the point about which all vehicle body force will be 
computed in later sections of this chapter. The convention used for the vehicle fixed 
frame, with unit vectors i, 7, K , is that its origin lies at the ship's center (origin at the half 
length), the x axis at main deck level, parallel to the longitudinal centerline, with the z axis 
vertically down. The vessel's center of gravity and center of buoyancy do not generally lie 
at the origin of the body fixed frame, nor are they collocated. These points are denoted by 
G and B respectively. The vector components of p, апа ñ, аге thus /x,/ - y, j +, J. 
and /[x,/+y,j/+z,k / . The location of the center of mass is important because 
Newton's Laws of Motion equate forces on a body to the rate of change of linear 
momentum ofthe center of mass and moments about the body center of mass to the rate 
of change of angular momentum. This is a particularly important point to bear in mind for 
ship and submarine motion dynamics as the center of buoyancy is determined by the shape 
of the submerged portion of the ship body while the center of gravity 1s determined by the 
distribution of the weight over the entire ship body. Having defined a coordinate system 
that will be used to describe a ship's position, there 1s a need to define angular orientation 
of the ship's body, leading to the definition of translational and rotational velocities and 
accelerations. 

3. Angular Position In The Global Reference Frame 

There are several different ways that the attitude of a vehicle can be described in 
reference to the global frame. The most usual and common method is to define three 
angles called Euler angles that uniquely define the angular orientation of the vehicle 


reference frame, relative to the global reference frame. One problem with the Euler angle 


approach is that a singularity exists when one of the angles reaches 90 degrees (an aircraft 
in pure vertical flight cannot distinguish its azimuth angle from its roll angle). This 
limitation which can sometimes, although rarely, cause trouble in flight simulations and 
control computations can be overcome by the use of quaternions which introduce four 
rather than three variables to describe angular position. In this presentation, however, we 
will use Euler angles as it is the most widely used method, although the use of quaternions 
has found favor in robotics applications and computer graphics. 

While any consistent definition of three base angles would be sufficient, the most 
convenient formulation of these angles is in the widely used military terms of an azimuth, 
elevation, and spin notation. The rates of change of these angles do not generally 
correspond to the other commonly used angular rates describing angular velocity of a 
vehicle, that is yaw rate, pitch rate and roll rate except. as will be seen, where motion is 
limited to small angle rotations. 

4. Rotational Transformations 

Vehicle attitude is important in defining position. Under dynamic conditions, the 
pitch or roll can cause problems for manned vessels and the heading is critical in 
navigation. For the purpose of considering angular rotations, consider a forward 
transformation from a coordinate triad aligned with the global reference frame and perform 
a sequence of three rotations to finally align the result with a frame that is assumed to be 
parallel to the vehicle body coordinate axes, and moving with the vehicle at all times. 
Begin by defining an azimuth rotation, y, as a positive rotation about the global Z axis. 
Next define a subsequent rotation 8, (positive up) about the new Y axis, followed by a 
positive rotation $, about the new X axis. The triple rotational transformation in terms of 
these three angles, is then sufficient to describe the angular orientation of the vehicle at any 


— 


tme. It follows that any position vector, К 


о? 


in an original reference frame given by 


К.-/х,.у,2,/, will have different coordinates in a rotated frame when an azimuth 
rotation by angle y, Is made about the global Z axis. 

If the new position is defined by, R, = /x,,),,z,/, it can be seen that there is a 
relation between the vector's coordinates in the new reference frame with those that it had 
in the old reference frame. Using trigonometrical relationships and Figure 2.2 as a guide, it 
follows that, 

X, = X cosy +Y, siny (2.1) 
y =— X. sing + Y coswv. (2.2) 
This relation can be expressed in matrix form by the rotation matrix operation, 

R z[1,;]R; (2.3) 
where the rotation matrix ШЕ represents an orthogonal transformation. 
Premultiplication of this rotation matrix with any vector, R,, will produce the components 
of the same vector in the rotated coordinate frame. Continuing with the series of rotations 
results in a combined total rotation transformation, 

(9,0, w) = T(9)1(0)T(w). (2.4) 


In expanded notation equation 2.4 takes the form; 


| cos y cos sin y cos — sin 
»" y sinOsinó—sinycoso sinysinOsind+cosycosd | cosO sino |, 


cosy sin6cosd+siny sind  siny sin8cosQ—cosy sino совӨсовф 
and it can be said that any position vector in an original reference frame may be 


expressed in a rotated frame with coordinates given by the operation, 


ког өк (2.5) 


Figure 2-2. Azimuth Rotation. 





y 


5. Kinematics 
Kinematics deals with the relationships of motion quantities regardless of the 
forces induced by their prescribed motions. Descriptions of a ship's position both 
translational and rotational, will need to be related to velocities, both translational and 
rotational, prior to extending the situation to accelerations. The connection between 
translational velocity and the rate of change of translational position is straight forward. 


Define the global velocity as, 


(2.6) 


Pel 
II 
Ку. не. ра. 


This vector will have components that are different when seen in a body fixed frame. 
Define the body fixed components of the global velocity vector as [u, v, w]. These 
components, in terms of the above global quantities are given by the forward 


transformation defined earlier to be, 


X 
v |- T(6.9,v) Y |. (2.7) 
у 2 


It is a simple reverse coordinate transformation from body fixed to global coordinates to 


see that, 


X и 
Y |=T'(0,6,w) v |. (2.8) 
Z и 


This shows that the progression of a vehicle in a global frame clearly depends on its local 


velocity components and its attitude. Put simply, u corresponding to a vehicle's forward 


speed (surge), v corresponding to a side slip velocity (sway), and w corresponding to any 
velocity component in the local Z direction (heave) and the vehicle's global velocity 
components depend on heading, pitch, and roll attitude. 

The connection between angular attitude and angular velocity is not quite so 
apparent. At first sight, it is tempting to define the instantaneous angular velocity of the 
vehicle simply as the rate of change of its angular position defined by the Euler angles. 
This 15 erroneous however, because the rotation 0, was defined as a rotation about the 
intermediate frame after a rotation y had been made. Vehicle inertial angular rates are 
defined in terms of components that have angular velocities about the global axes. It is 
necessary to relate both Euler angles and their rates of change to angular velocity 
components about the global axes to their components lying along the body fixed axes in 
any attitude. The prime reason for this is that it is difficult to construct physical sensors to 
measure rates of change of Euler angles. However, rate gyros in common use today are 
easily constructed to measure the components of the inertial angular velocity of a vehicle 
that lie along the vehicle's body axes. It follows that the instantaneous angular velocity of 
the vehicle can be related to the instantaneous rate of change of angular orientation only 


after considerations of the intermediate transformations used. In other words, if œ is 


defined as the angular attitude vector, « =[w,0,o]; and the inertial components of the 
vehicle angular rate lying along the body axes as @ = ра then, & = f (w). The details 
of the nonlinear functional relations involved are provided by viewing the rate of change of 
the rotation V as a vector quantity lving along the original Z axis. The rate of change of 


the angle 0 is viewed as a vector quantity lying along the Y axis of the first intermediate 


frame, and the rate of change of the angle ф 15 viewed as a vector lying along the X axis of 
the final (body fixed) frame. Each of these component rates of change of angular position 


has component parts that project onto the final frame and it is the sum total of all the 
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components that give the total angular velocity as seen in the final frame of reference 


Using the required transformations for the rate components from each Euler angle, we get, 


p 0 0 ф 
q |= T(#)T(0)T(v) 0 |+T/($)T(9) 0 |+ T() | 0 |, (2.9) 
r V 0 0 
with the result; 
p -y sin0 +6 
qi- ц/сіпӨ--Өсозф |. (2.10) 
r ў со50со5ф – Ө ѕіпф 


Inverting equation 2.10 yields a solution for the rates of change of the Euler angles in 


terms of the body fixed components of the angular velocity vector, 


ф p+qsingtan8+rcosotan® 
6 |- qcosQ —r sino | (2.11) 
V (азіпф--ғсовф) / соѕӨ 


Notice that for small angular rotations, as expected, 


ф 
8 
у 


І. 
асы 


At this point the kinematic relationships between velocities, as seen in the body fixed 
frame, and the rates of change of global positions and Euler angles have been defined. The 
resulting set of differential equations forms a consistent set in that given a set of vehicle 


velocity data versus time, its position and attitude may be computed. 


6. Translational Equations Of Motion 
The global acceleration of the center of mass is derived by differentiating the 
velocity vector, В, taking into account that the center of mass lies in a rotating reference 
frame. Considering the total differentiation, the global acceleration of the center of mass 
becomes, 

В, =” +@хр +@хф@хрб„+ф@х, (2.12) 
where v К. The translational equation of motion is found by equating this acceleration 
to the net sum of all forces acting on the vehicle in three degrees of freedom (X,Y,Z). One 
important factor to recognize is that the equation of motion derived in this manner is a 
vector equation with the components expressed in the body fixed frame and unit vectors 
i, jand k. This has been deliberately done because the dominant forces acting on a 
submerged body in motion are developed in relation to the shape of the vehicle and are 
more conveniently expressed in relation to the body axes. Equating the applied force 
vector to the acceleration, results in, 

F =m{v+0xp,+@x@xp,+@x,}. (2:15 
The applied force vector is composed of gravitational (weight) and hydrostatic (buoyancy) 
forces together with any hydrodynamic forces arising from relative motion between the 
vehicle and the ocean water particles. These will be discussed in further detail in following 
sections 

7. Rotational Equations Of Motion 
To develop the rotational equations of motion, the sum of applied moments 
about the vehicle's center of mass is equated to the rate of change of angular momentum of 
the vehicle about its center of mass. This equating will provide the necessary equations of 
motion for the remaining three degrees of freedom. In the practical case of marine 


vehicles, however, the statement just made is modified slightly because it 1s much more 


[2 


difficult to assess the vehicle's mass moments of inertia about its center of gravity (CG) 


As the CG changes with loading; it becomes simpler to evaluate the mass moments of 


inertia about the body fixed frame which lies along the axes of symmetry ofthe vehicle in 


most cases. The mass moments of inertia for the vehicle to be computed are; 


i ху x 
а ы) 
EX 
The angular momentum of the body is thus; 
DAE 


resulting in the total applied moments about the ongin given by, 
M, = iun tpa x (mh ). 


Since the rate of change of angular momentum is given by, 


Ф 
==» 


Н„=1®+@хН,„, 


and the acceleration of the global position vector is given by, 


— 


R =0+0xi, 


о 


then the rotational equation of motion in vector terms is given by, 


М, = 1@+@х(1,®)+т{р„ ху +р„ х®х}. 


(2.14) 


(15) 


(2.16) 


(2.17) 


(218) 


(2.19) 


The applied moments about the body fixed frame arise from a static balance of 


weight and buoyancy effects to achieve the proper trim and heel, and hydrodynamic 


moments from the forces applied through appendages such as control fins, hydrodynamic 


effects from waves, and hydrodynamic effects from relative motion between the vehicle 


and the water. 


At this point, there are three translational equations obtained from equation 


2.13, three rotational equations obtained from equation 2.19, and six unknown velocities 


(Ф and v). In itself this is a consistent set of dynamic equations if the weight and buoyancy 
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terms are always self canceling. However, with weight and buoyancy being applied in a 
global vertical direction, and also being applied at different locations, they represent forces 
that are dependent on the attitude of the vehicle. Notice that without weight and buoyancy 
forces, the six consistent equations expressed in the body fixed frame of the vehicle could 
be solved for the vehicle's velocity (given applied force descriptions). With weight and 
buoyancy acting, we need additional equations (constraints) that will link vehicle attitude 
to motion so that the combined set of equations can be solvable. The constraints to be used 
are developed using the Euler angles which are a general set of angles used to define the 
attitude of the mgid body. From these angles a relationship between rate of change of 
attitude and the body rotational rates given earlier as the angular velocity vector 
Q — / p.q.r J, can be developed. 
8. Incorporation Of Vertical Forces Into The Equation Of Motion 

The weight and buoyant forces that act at the centers of gravity and buoyancy 
must be defined from static analyses. For submerged bodies the weight and buoyancy force 
vectors do not change with vehicle attitude. For a surface ship, B will change with 
attitude, and the mghting moments must be computed from Naval Architectural 
considerations. Assuming that weight and buoyancy are fixed in relation to the body fixed 
frame, W =0/ +0/+(mg)K, and B=0/ +0/-(pV)K. Since the weight and buoyancy 
terms in the applied forces act in the global vertical direction, they must be transformed 
into components in the vehicle fixed frame before they can be added into the equations of 
motion Returning to equation 2.4, we therefore find the components along the vehicle 
fixed frame to be the third column of that transformation matrix. The net vertical force 


components become, 
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—sin® 


7, З (Й - ВД созд5інф |. (2.20) 


сов Өсовф 
The weight portion of the vertical force acts at the center of gravity of the vehicle. The 
buoyancy portion of the vertical force acts at the center of buoyancy. Because these forces 


act at positions away from the body center they exhibit a moment about the body center 


given by, 
— 5110 — sin 0 
m, =Wp, x| cosOsing |- Bp; x} cosO sing |. (2.21) 
соѕдсоѕф cos0cos 


This moment will not be zero even if W and B are identical because it is not usually the 


case that Б. and p, are collocated. For static stability it is advisable to locate the center of 


gravity below the center of buoyancy The total vertical force vector can be written as, 


а 


and is added negatively to the left hand side of the equations of motion. 
9. Development Of The Full Six Degrees Of Freedom Non-Linear Equations 
Of Motion For A Marine Vehicle 

Define a vector X of vehicle body frame velocities to be, ¥=/u,v,w, p,g,r/, 
and a vector Z of global positions to be 7=/X,Y,Z,0,0,w/, then considering M as a 
6x6 mass matrix including translational and rotational inertial elements, the equations of 
motion can be written in the following vector form, 

My + f(x)+ F(z) = E, (2.23) 

and, 


7+ e(x,2) =0. (2.24) 


[5 


Therefore with suitable knowledge of the excitation force and moment loads, 


F, as a function of time and/or vessel motion, a solution for the vehicle's dynamics can be 
obtained. A more detailed insight into the development of the twelve differential equations, 


in first order form given by the foregoing analysis shows, 


X 


/ 
т{ў + © х ре} + т(®х ® хр +@х)+ }(:)=| Ү, |, (2.25) 
2, 
апа, 
К, 
Là тр. Х ў) ФХ (16) тр. Хбх ў й, (2) а | М, |. (2.26) 
М, 
It helps here to define the cross product coefficient matrix so that, 
(0x, 2 -p, xà - -[cros(p,) à (2.27) 
where, 
0 d EE D 
СТОР 79/5 [. (2.28) 
OR eG 9 


Now collecting the mass matrix coefficients into a 6x6 matrix including the inertia cross 


coupling effects, results in the following equation; 
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т O O му = 
0 m 0 т- 9 хы 
e 0 0 m D eeu as (2.29) 
0 24, e CO - 
| E 0 -х 1 x D т | 
= 0 ВИ | 


The remaining terms on the left hand side of the equations of motion arising from the 


centripetal and coriolis accelerations become, 


m(QxQ xps * Oxv) | om 


по) rea ali x Q) x v ) 

The double vector cross products are nonlinear in the primary velocity variables and hence 
the need for the nonlinear functional, /(е) The reader can perform the indicated 
manipulations to express individual equations within the set if so desired. 

The components of the hydrodynamic and external forces and moments acting 
on the vehicle body are separated in the above analysis into six components each acting 
along the vehicle body fixed coordinate axes and form the total vector of forces and 
moments as, 

Et) m [Xo(t) Y) ZR) M, (NS, (2.31) 
where the vector components in order refer to the surge, sway, heave forces, and the roll, 
pitch, and yaw moments respectively. 

The long form of the six degrees of freedom equations of motion can be written 
as follows, 


m[u-vr *wq -x,(q^ ^r )*t ys( pq -r)* z;, (pre q)] - 


(2.32) 
оо. 
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m[ v -- ur - wp * x«( pq er)—- ys( p. *r )*zs(qr-p)] - 


2.33 
(Й – В)соѕ0 ѕіпф + Р, ( 


m[w-uq-*vp*x4( pr -q)* ys(qr t p)-z,(p +g )]= 


| (2.34) 
(W — B)cos8coso * Z, 


|р -L)gr-L(pr-d)-1,.(q -r)-lIj(pq*r)^ 
m[ y,(w —uq t vp)- (Ў t ur -wp)] - ( yaW — y,B)cos0 coso (2:55) 
-(zgW —zg B)cos0 sino - K, 


[q+(1 =l ə)pr= l (агер) (раг) 1 (ранее 
m[ x,(w -—uq +vp)—z (ü—vr+wq )]J=-(x.W -х,В)с050с05Ф (2.56) 
-( z;sW — z, B) sin0 * M, 


Lp дра Е ОЛА ага 
m[ x.(»- ur -wp)— ys(u — vr -*wp)] Z( x;SW - x,B)cos0 sin (2287) 
*( y9W — yg B)sin0 * М, 
while the kinematic relations for the vehicle rate of change of attitude and motion over the 
ocean bottom require equations 2.11 and 2.8. 
10. Adaptation Of The Non-Linear Equations Of Motion To The Vertical 
Plane 
Restricting the motions of the vehicle to the vertical (dive) plane, the only 
significant motions that must be incorporated to effectively model the vehicle in the dive 
plane are, the surge velocity (и), the heave velocity (w), the pitch velocity (q), the pitch 
angle (O) and the global depth position (Z). This restriction simplifies the twelve 
previously developed equations to a system of four non-linear equations of motion, which 
are; 


т(№ ид год ход)“ 7, (2.38) 
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l,d* mz,wq - mx,(w -ug)- M, 
8-4 
7 = -usin0 +w cos9 


where, 


2,-2,4%2,М%2,4 %2,им- 


>р | (Бері E HD aes Qr - B)cos@+u° (2,0,%2; б,) 


tail | | 


and, 


M .=M.qg+M.w+M ца М им 
фнди 
1 e Cbl xq) хах 


tail 


š 2 
- (xW -xgB)cos8 - (z-W - z ,B)sin 0+u (M, ó, + М,, 6, ) 


(2.39) 
(2 40) 


(2.41) 


аказ) 


results from expanding the Z, and the M, terms from (2.31) in a first order Taylor series 


expansion and incorporating both hydrostatic and fluid drag forces. 


Equations 2.38, 2.39, 2.40 and 2.41 can be linearized for a level flight path when 


the dive plane angle is zero, ie. 0, 20. By setting all time derivatives to zero, and 


neglecting for the moment the hydrodynamic drag terms the following are obtained; 


Z.uw * (W — B)cos0-0 
Ва Ни ух 2D) COSO=( 20 =o, D) SIG = 0 
q= 0 


—usinO 4 исоѕ = 0. 
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(2.44) 
(2.45) 
(2.46) 


(2.47) 


If the assumption is made that the vehicle is neutrally buoyant, then it can be said that 
x, =x, and W - B. From this equations 2.44 and 2.45 can be reduced to; 

Z uw =0 (2.48) 

MM —(z ару ВО 0, (2.49) 

which yield the nominal position w, — q, 2 0 and sin6, — 0, resulting in the 0, solution as 


either 0, =0 or 8, = x. Choosing to linearize around the nominal point 6, = 0 results in; 


4 -244-0 (2.50) 
масмағамт0 (2.51) 
5тӨ = соѕ0 Ө = Ө (2.52) 
wcos0 2 (—w,sin0,)0-(cos0,)w 2 w. (2859) 


The linear equations of motion are then written as; 
(m— Z,)w -(Z,* mx; )q — Z,uw t (Z, * m)uq * 2и% (2.54) 
(1. — M,)d - ( M, * mx, )w Б М ше c (M, — mx, Jug - (z; - z, WO Mud (2:55) 
Ө-4 (2.56) 
2--ибем. (2.57) 
As equations 2.42 and 2.43 show, both Zs and Ms are a linear combination of 
the respective stern and bow hydrodynamic control surface coefficients and the respective 
input value of д. This makes the system of equations as a multiple input system. To 
reduce this system into a single input system the linear combination of control inputs will 

be modified into the following form; 


Z- (2, +02, ). (2.58) 


This will allow a single input 6 to control both stern planes and bow planes, and will 
cause the bow planes to be slaved to the stern planes. This technique is known as dual 
control. The value € will range from -1 to 1. The selection of the value of O. will allow 


the planes to operate as desired for the particular maneuvering condition, 1.e., & = 0 for no 
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bow plane control, œ = -1 for bow plane and stern plane control opposed to each other, 


yielding the maximum pitch moment, and & — 1 for bow and stern plane control in the 


same direction, yielding the maximum heave force. 


In state space form these equations can be written as; 


0 0 0 1 09 0 
і a-z au au Ollw bu 
w f: 1 12 + l ‚|6. 
4 а ана иона |b 
7 —H ] 0 0127 0 


where the coefficients a, and b, are given by; 
D, Z(m- Z,)( T, - M.) - ( mx, * Z, )( mx, * M.) 
a,D, 2 (1, - M,)Z, * (mx, * Z,)M, 
a,,D, =(1,-M, )(m+ Z,)+ (mx, + Z,)( M, —mxç ) 
Gu =-(mx,+Z,)W 
Бр, = (T, - M, )Z, + (mxç + Z, )M, 
a, D. = (m- Z. )M. + (mx + M. )Z_ 
a,,D, 2 (m- Z,)( M, -mxç)+(mxç+ M, )(m+ Z) 
a4,D,--(m-Z,)W 
раза ўт Z, )M, + (mx. + M. )Z, 


and z;4 2 z; — z, is the metacentric height. 


(2.59) 


(2.60) 
(2.61) 
(2.62) 
(2.63) 
(2.64) 
(2.65) 
(2.66) 
(2.67) 
(2.68) 


This linear system of equations becomes the basis for the control law design in 


the following section. 


B. CONTROL LAW DESIGN 


1. Introduction 


The control design problem can be stated as follows: Given the system; 
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x 2 Ax 4 B6, (2.69) 


where the state vector equation 1s, 


x = (2.70) 


0 
М 
a 
Z 
how do we find ò, such that the system will behave as desired. The type of control that is 
of interest in this problem is closed loop control, where ò is a function of the state x. 
Since the state x is used to determine the control effort ó( x ) it is called feedback control. 
2. Pole Placement 
A linear full state feedback control law is introduced in the form 
б=-Ах, (2.71) 

where K is the feedback gain vector to be determined such that the closed loop system of 
equations 2.69 and 2.71 has the desired system dynamics. Substituting equation 2.71 into 
equation 2.69 yields, 

xz(A-BK)x. (2 72) 
The actual characteristic equation of this closed loop system is given by 

det|A - BK — sI ] - 0. (2.73) 
The gain vector K can be chosen such that the actual characteristic equation assumes any 
desired set of eigenvalues. If the desired locations of the closed loop poles are chosen at 

= s for i-l,...,7, the desired characteristic equation becomes 

(5 — 5, )(5 — 5. ss О (2.74) 

The required values of K are obtained by matching coefficients in the two polynomials of 


the actual and desired characteristic equations. 
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Consider the previously developed linear state matrix equation 2.59. 


Substituting equation 2.71 into equation 2.59 results in the closed loop system 


4 Š 0 1 0 Ө 

x _| азо -hwk a-bwk a-bwk, бик |и (2.75) 
; Ааа ра, ася - bark, а.и- Б.Е, Бик, q | | 

7 =u | i ў < 


The characteristic equation of the closed loop system is given by 


_ 0 l 0 
2 2 2 
E азва б а НЕА Е a ubu ka Sbu ki T (2.76) 
2 2 2 жог; = 
алга - А a,,u-bauk, au-buk-s Ыш, 
-и 1 0 eS 


which after some algebra reduces to 


ая Вов Eds 


+(-С С, - C,k, - E )s + ( CDk, - Dk,) 20 e 
where , 
4 = B. =b (2.78) 
А, =-В = Би" (2.79) 
B, = (а,Ь -а.Ь, )и (2.80) 
В (а раи (2.81) 
Сер (2.82) 
C, 2 (a,b, * b, - a,b, )u' (2.83) 
D, = (a, b, — a,b )u* (2.84) 
E, = (a qu (2185) 
E gis adsensu: (2.86) 
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B Tq аа ие (2.87) 
Now if the closed loop poles are placed at —p,, — p,, – р,, апі – р,, then the 


desired characteristic equation is, 


(STD) St pols үре (2.88) 
ОГ 

5% +05 0s casa, = 0 (2.89) 
with, 

OQ; -Dtptp,tp, (2.90) 

©, = рр; + РР; + РР. + Р.Р, + Р.Р. + Р.Р. (2.91) 

0; = Р.Р; + РР:Р. + РРР, + PPP (2.92) 

=, = рр. РР. (2.93) 


The control gains can now be computed by equating coefficients of the actual and desired 


characteristic equations, 


Ak, + Ak, 2 -a, - E, (2.94) 
pk + 5.k EE —O0 CE. (2:95) 
G,k, + C,k, + C,k, = w, + Ë, (2.96) 
(Di + Do )k4 = Q4. (2.97) 


Now that a method for computing the gains of a controllable single input system 
that will place the poles at any desired location has been developed, the selection of pole 
location must be accomplished. 

3. Pole Location Selection 

In a typical second order system control law design the transient response 
specifications will be given. This results in an allowable region in the s-plane where the 
desired location of the poles can be obtained. For higher order systems the concept of 


dominant roots can be employed. In selecting poles for a physical system it is necessary to 
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look at the physics of the system. If the poles are specified too negative, a very small time 
constant for the control system will result, and the physical system may not be able to 
react that fast. The control law 4 2 —Kx, implies that for a given state x the larger the 
gain, the larger the required control input. In practice there are limits typically placed on v 
(actuator size and saturation). Occasional control saturation is not serious and may even 
be desired. A system that never saturates is in all likelihood over designed. 

Considering the control law design to stabilize the submarine to a level flight 
path at 0 20 it will be required that the submarine return to level flight, after a small 
disturbance in 0 or Z, within the time it takes for the vehicle to travel three ship lengths. 
Since the submarine is about 14 feet long and it travels at 5 ft/sec, the required recovery 
time is about 10 seconds. This means that the time constant is about 3 seconds, and the 
closed loop poles should be placed at approximately —0.3. 

Control design using pole placement is very easy using MATLAB. The 
appropriate MATLAB command is place, which accepts as inputs the A and B matrices 
along with a vector of the desired closed scis poles, and returns the gain vector K. Using 
equation 2.59 with a nominal speed of 5 ft/sec and the vector 
p-[-9.3 -0.31 -0.32 -0.33/, (place does not like poles in the exact same location), 
the gain vector K was calculated using MATLAB and by simultaneously solving equations 
2.94 through 2.97 yielding the same results. Substituting the gain vector and state variable 
vector into equation 2.71 results in a control law of 

ò = —(—0.99170 —0.8333w — 0.6026q 4 0.0351Z . (2.98) 

Using this control law along with equations 2.38 through 2.43 a simulation 

program was developed to investigate the vehicle's response to initial disturbances, the 


results of which will be presented in the following chapter. 
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HI. PROBLEM IDENTIFICATION 


A. VEHICLE SIMULATIONS 
With the control law developed and the vertical plane equations of motion defined, a 

vehicle simulation program was developed (see Appendix A). This program was used to 
simulate the vehicle's response at a variety of speeds and initial disturbances. The 
simulations were conducted with several simplifications applied to the vertical plane 
equations of motion (E.O.M.), equations 2.38 through 2.43. The simplifications or 
constraints applied to the E.O.M. were ; 

1) to consider the body drag forces negligible, 

2) to consider the vehicle to be neutrally buoyant, i.e., (W=B), 

3) to assume that the locations of the center of gravity and center of buoyancy, 

in the X-Y plane, are coincidental, and 

4) that the rudder is restricted to + 23 degrees. 

These assumption resulted in the reduction of equations 2.38 through 2.43 to the 


following system of equations; 


l 0 0 Olé q 

DU mE 2А Ol w т(ид + год" )* Z,uq * Zaw + u` Z, (3.1) 
о =M, и i m( -z,wq )* M, uq * M,uw — Bigg sin® + wM; 

0 0 0 117 —u sin0 + w cos 9 


forming the basis of the initial vehicle simulations. 
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The vehicle's response to small disturbances, as stated earlier, was investigated over a 
wide range of speeds. A small sample of these simulations is shown in Figures 3-1 through 
3-3. 

As Figure 3-1 depicts, the vehicle returns to level flight in approximately 10 seconds. 
This response is as expected based on the control system design of Chapter IL. The 
simulation results represented in Figure 3-2 also show the vehicle steadying out in a level 
flight at the desired depth, however the vehicle requires a significant amount of time for 
this to occur. The reason for this extended amount of time is that the control law gains 
were set for the nominal operating speed of 5 fps, and the vehicle is being simulated at 
one-third of that speed. This results in a very slow return to ordered position caused by 
the reduced hydrodynamic effects at this speed. 

The final plot of the sample vehicle simulations, Figure 3-3, illustrates a problem with 
the vehicle stability. As can be seen, the vehicle does not return to the desired position, but 
instead steadies out at a pitch angle of eight degrees, with a steady state dive plane value 
of 23 degrees. This same type of E with different pitch angle values, was 
observed in all the simulations conducted below approximately 1.9 fps. These results 
indicate that there is a critical point, or speed at which the vehicle stability is affected. 


Determining this point will be the basis for the remainder of this chapter. 


B. CRITICAL SPEED IDENTIFICATION 
1. Eigenvalue Analysis 
Consider the nonlinear system of state equations 3.1, which can be written in the 
form, 
ил (3.2) 


It is known that the equilibrium points, x,, of the system are defined by 
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Figure 3-1. Vehicle Response For a Nominal Operating Speed of 5 fps. 
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Figure 3-2. Vehicle Response at 2.0 fps. 
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Figure 3-3. Vehicle Response at 1.885 fps. 
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f(x,)=0. (3.3) 
This is a nonlinear system of algebraic equations and it may have multiple solutions in x,, 
which means that the nonlinear system of equations may have more than one position of 
static equilibrium. If one equilibrium value x,, is selected, the stability properties can be 
established by linearization. The linearization converts equation 3.2 into the form 
Х- Ах, (3.4) 
where A is the Jacobian matrix of f/x) evaluated at the equilibrium point x,,, 


г 
= р (3.5) 


Хо 


and the state x has been redefined to designate small deviations from the equilibrium x,, 

x— x-—x.. (3.6) 
As long as all eigenvalues of A have negative real parts, the linear system will be stable. 
This means that the equilibrium x, will be stable for the nonlinear system as well. This 
statement is nothing more than Lyapunov's linearization technique. 

The question that must be answered is, what happens if one real eigenvalue of 
the linearized A matrix is zero? The interesting case here 1s when the rest of the 
eigenvalues have all negative real parts, otherwise x, is unstable and the problem is solved. 
If the case of a zero eigenvalue appears too specialized to be of any practical use, consider 
this: Assume that f(x) depends on one physical parameter and that physical parameter is 
allowed to vary over some range. Then clearly A will depend on that parameter and as the 
parameter varies, it is possible that one real eigenvalue of A will become zero for a specific 
value of the parameter. The problem is then to establish the dynamics of the nonlinear 
system as one real eigenvalue of A crosses zero, i.e., goes from negative to positive. As 
the solution evolves in time, things are interesting only along the direction of the 


eigenvector that corresponds to the critical eigenvalue. Along the rest of the directions in 


3l 


the state space, everything should converge back to the equilibrium; remember it was 
assumed that all remaining eigenvalues of A have real negative parts. Although, strictly 
speaking, this is a true statement for linear systems, there are technical reasons that force it 
to be true for nonlinear systems as well, the only difference is that the corresponding 
directions in the state space are curved instead of straight. 

By taking the previously developed characteristic equation 2.77, and 
remembering that the roots of this equation are the eigenvalues of interest, it is easy to see 
that if ©, of equation 2.93, is set equal to zero, then one eigenvalue will be zero. If this is 
done, and recalling that & , holds a non-zero value. equation 2.97 reduces to 

(a b, -a,b и? +2.„(а.Ь, -а,„Ь ) = 0. (3.7) 

Since all parameters of equation 3.7 have a fixed value with the exception of v, there must 
be some value of и that causes the linear A matrix to become unstable. Recognizing this 
fact will allow equation 3.7 to be solved in terms of u yielding 

и = E | (3.8) 

(a,b, —a, b) 
Using equations 2.60 through 2.68, equation 3.8 can be simplified into the following form; 
1/2 

Ж bus E | nd 
Substituting the appropriate values from Appendix B, equation 3.9 can be solved for the 
value of u that causes the systems to become unstable. This results in a value of u = 
1.8979 fps using an w = 0. This value corresponds very closely to the value of 1.9 fps 


obtained by vehicle simulations. 


As a check to this solution, a MATLAB program was written to compute the 
closed loop eigenvalues of the matrix equation 2.59 as a function of the speed и. 
Figure 3-4 shows a plot of the real parts of the eigenvalues versus speed. It can be seen 
that the real part of the eigenvalue, represented by the dash-dot line, becomes positive at 
speeds less than approximately 1.9 fps supporting the earlier findings. (For clarity 
purposes, the values of the eigenvalue represented by the dash-dot line were scaled by a 
factor of ten.) 

2. Steady State Analysis 

To determine the existence of multiple steady state solutions, for the system 
represented by equation 2.59, it 1s necessary to perform a steady state analysis on the 
system of equations which models the vehicle dynamics. Using matrix equation 3.1, and 


setting all time derivatives equal to zero, results in the following set of equations; 


4-0, (3.10) 

Z uq + Z ¿uw +u°Z,ó = m(—uq -zgq' ), (3.11) 
Мид + M.,uw — Bzos sin9 + 1 M,Ó = mz wq , (3.12) 
-usinO t wcosO - 0. (3.13) 


Simplifying these four equations yields the following two equations in terms of 0, и, апа 
the physical parameters of the vehicle; 
Z.u' tanü* w Z,620 (3.14) 
M wu’ tan® - Bz, sin6+u°M,6 = 0. (3.15) 
Now if equation 3.14 is multiplied by M,, equation 3.15 1s multiplied by 2, and the two 
resulting equations are set equal to each other and simplified, the following single equation 


Is obtained; 


[(Z, M, - M,Z,)w * Z,Bz,, cos®]sin® = 0. (3.16) 
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Figure 3-4. Real Part of Closed Loop System Eigenvalues as a Function of Speed. 
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Through inspection it can be seen that equation 3.16 could have multiple 
solutions in Ө, besides the trivial solution of 8 = 0. So, if this equation is rearranged and 


solved for 8 + 0, solutions occur when 


| w(M,Z,- Z, M,) 


050 (3.17) 
25В:св 
Therefore, the steady state value of O is represented by 
0. cor MLZ -E My d (3.18) 
Zo 


which may have multiple solutions based on the value of the speed uw. As discussed 
previously, these multiple solutions were evident in the simulations conducted earlier in 
this chapter. 

Knowing that the maximum value for cos0 is equal to one, the right-hand side 
of equation 3.17 will be true for values less than or equal to one. If this constraint 15 


imposed on 3.17, then the equation may be manipulated into the following form 


р ee уе (3.19) 
(M,Z, ш Z,M;) 


Rearranging this equation and solving for the upper limiting case yields the same results as 
equation 3.9. Again this supports the critical speed values discussed earlier. 

Equation 3.19 can be expressed in a non-dimensional Froude number form by 
dividing the equation by both the gravitational constant g and the physical parameter zog 
and taking the square root of the result. This will convert equation 3.19 into the non- 


dimensional form given below; 


Pre | ШЕ se кой э (3.20) 
87GB 8(М,%-2,М;) 








With the steady state O solution derived, equations for both the steady state 
value of 5 and steady state value of Z can be determined. Using equation 3.14 and 
substituting the steady state value of 0 will allow the steady state Ó equation to be 


obtained; 





о. = = tan... (3.21) 
Z 


5 
To achieve the steady state solution of Z, begin by writing equation 2.71 in 
expanded general form as 
6=—-(kO+kw+kgqtk,Z). (3.22) 
Then by applying the same conditions to equation 3 22 that resulted in equations 3.14 and 


3 15, and using the previously developed steady state equations for O and ò, an equation 


of the form 





2, апе, +0, + А.иіапӨ,, 
2. = ст ss (3352.3) 
15 obtained. 

The loss of stability of an equilibnum and the generation of additional 
equilibrium states, is called a pitchfork bifurcation and is very common in nature. The 
buckling of a beam is one such example. Using equations 3.18, 3.21, and 3.23 as a basis 
for another MATLAB program, the steady state solutions of 0, ô, and Z versus Froude 
number were investigated with the results displayed in Figures 3-5 through 3-7. These 
three plots are referred to as supercritical pitchforks, so named because upon the loss of 
stability of the trivial equilibrium the additional nearby equilibrium states are stable. 
Graphically, this can be represented in Figures 3-5 through 3-7, where the solid curves 


represent stable and dotted curves represent unstable equilibria [Ref. 17]. These three 


figures also demonstrate the control input saturation. Upon control surface saturation, the 
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steady state equations for 0, ó and Z are no longer valid due to the limit placed on the 
maximum plane angle. The steady state equations that hold for this region of saturation 


are given below; 


6, = + 22.9' (0.4 radians ), (3.24) 
мош--- 3 25 
55 2. ( ) 
and, 
0. = «ін! Mau ð, – М,им,, | (3.26) 
28 


As can be seen in Figures 3-5 and 3-6, at an approximate Froude number of 1.05 the 
control surfaces saturate resulting in the linear solutions displayed below that saturation 
Froude number. In Figure 3-7 there is no steady state solution because if the steady state 
values of 5 and @ are placed into the Z equation, Z is non-zero below the saturation 
Froude number. 

A parameter introduced in Chapter II was X, which allowed us to go from a 
multiple input system to a single input system. This parameter has a significant effect on 
the location of the critical speed and/or Froude number. Figure 3-8 shows the relationship 
been the critical speed u,, and the metacentric height zog for different values of a. 

These three curves shown in Figure 3-8 can be converted into a single curve by 
plotting & versus Froude number, апа is shown in Figure 3-9. As can be seen in this plot 
the critical Froude number varies from 0.81 to 1.22 over the allowable range of C. 

In addition to the effect that & has on the location of the critical Froude number, 


it also effects the magnitude of the steady state values at a given Froude number and the 
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Figure 3-5. Steady State Values of 0 vs Froude Number. 
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location at which the control surfaces saturate. When comparisons are made, it can be 
shown that as the value of © becomes more positive the location of the critical Froude 
number and steady state values of 0 and Z become greater. Comparing the changes in the 
ò steady state solutions, the shape of the ô pitchfork does not change however, the 
location of the critical point moves in the direction of increasing Froude number as С 
increases. The results of this discussion are shown graphically in Figures 3-10 through 
3-12. 

An interesting result is displayed in Figure 3-13. This figure depicts the 
relationship between & and the critical or saturation Froude number. At the lower Froude 
numbers, there is very little difference between the critical Froude number and the 
saturation Froude number. This demonstrates that upon reaching a low critical Froude 
number the control surface immediately saturates attempting to keep the vehicle stable. At 
the higher Froude numbers there is a significant difference between the critical Froude 
number and the saturation Froude number. This is because as the Froude number 
approaches 2.78, the value used in the control law design conducted in Chapter II, the 
control surfaces will not have to saturate to maintain stability. 

3. Controllability Analysis 


In general, every system of the form, 


х = Ах+ Ви, 
у= Сх, 


(5227) 
can be divided through a series of transformations into four subsystems: 

1) A controllable and observable part. 

2) An uncontrollable and observable part. 


3) A controllable and unobservable part. 


| 4)  Anuncontrollable and unobservable part. 
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Figure 3-10. Comparison of Steady State Values of 0 for Different Values of c. 
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This is known as Kalman's decomposition theorem. Now, since the transfer function of 
any system is determined by the controllable and observable subsystems only, then the 
transfer function may contain less information than is actually needed to model the 
complete system. The precise definition of controllability 1s: 
A system is said to be controllable if any initial state x(t,) can be driven to any 
final state х(1,) using possibly unbounded control u(t) in finite time 1,< f < ty 
From the state equations 3.27, this should depend only on Æ and B, where A isan X n 
square matrix. The criterion for controllability is as follows: Compute the controllability 
matrix 
C=) BAB AB в. (3.28) 
and the system is controllable if and only if the rank of c (the number of linearly 
independent rows or columns) is 7. Roughly speaking, c shows how possible it is to 
change the state of a system using the input. For a single input system B is n X l and c is a 
square matrix. The test is then that c be nonsingular 
det(c) #0. (3.29) 
Now, if the controllability matrix is formed using the A and B matrices from 


equation 2.59 then the following matrix results; 


MU. x (3.30) 


С Ср Ca Си 


where, 
C= О, (3.31) 
с. = Би, GED) 
e a Dea ape (3.33) 
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араасаа а и) рика с заа иаи). (3.34) 


с, ари”, (2 39) 
_ 3 3 
Gd ub tanu Db, _ (3.36) 
ери ау и ва ам Дм (ас, аа” ta ан), (3.37) 


cup: (аа) u` +a a u )+a,u(a,z ов + ана, зи" +a a,u” ))+ 
Б.и? (а, ова + а 432и +а,и(ар и ааш 1) (3.38) 
а иа ге каа ци“ нати“ )), 


Grabs. (3.39) 
Ca u b sn pb (3.40) 
с, = Би (a, a. aa, уйи (асылын ры ну (3.41) 


Dur (a u(a a u" taca ii Pa ER м? ))+ 


bau? (qr. gd а ага а иѓа а и" рага”) (3.42) 
а (аг, ага” За), 
сі 50, (3.43) 
жылы (3.44) 
c, 2 au b, t (a.u -u)b,w, (3.45) 
and с. = Би? (а, ди? + a.u(a,u-u))+b.u`(a,z +a a.u`+a,u(a,.u-u)) (3.46) 


Now, if the determinant of the controllability matrix is computed as a function of 
Froude number, then it can be shown that the system becomes uncontrollable at the critical 
Froude number. The graphical results of this computation are shown in Figure 3-14 for 
different values of the parameter Q. 

With the critical parameter that causes the vehicle to lose stability identified, 1.e., 
the Froude number, and the effect that the parameter o has on the location of this critical 
value understood, the simplifications discussed at the beginning of this chapter can be 
removed from the E.O.M. and an analysis on the effects of incorporating these additional 


parameters can be conducted. This detailed analysis will be the topic of Chapter IV. 
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IV. BIFURCATION ANALYSIS 


A. ASYMMETRIC PITCHFORK 

Now that the critical parameter, the Froude number, has been identified, the 
simplifications applied to the E.O.M. in the previous chapter can be removed. The 
resulting system of equations is given by equations 2.38 through 2.43. To determine the 
existence of multiple steady state solutions for this new system of equations it is, once 
again, necessary to perform a steady state analysis as conducted in Section B.2 of the 


previous chapter. This analysis will reduce the system of equations into the following; 





Z, uw - -pC; Aw (Й - B )cos8 + 2,и°ӧ = 0 (4.1) 


Мим = PCox Aww - (XW - x,B)cos0 - ( z,.W —z,B)sin0- M,w 8 = 0, (4.2) 
with equation 3.13 also obtained from the steady state Z equation. By multiplying equation 
4.] by M; and equation 4.2 by Z, and setting the resulting equations equal to each other 
the following equation is obtained; 

] 
(Z,M,- M,Z, jw - (7 pAC, JC My — x ,Z, )wiw|* (43) 
((W—B)M, *( x,W —x,B)Z, )cos0 * Z,(z,W — z, B )sin@ = 0 
This equation, along with, 
w -utan8, (454) 


from equation 3.13 is the exact steady state solution set with all parameters included. 


5] 


Figure 4-1 shows the results obtained when equation 4.3 is solved numerically, using the 
FORTRAN program in Appendix A, for a small value of x,, (the difference in the location 
of the center of gravity and buoyancy), with the submarine neutrally buoyant (W = B). 
This small value of x, perturbs the pitchfork, as seen, from the symmetric pitchfork 
displayed earlier in Figure 3-5. This comparison shows that the critical or bifurcation point 
has moved to a lower Froude number, with the stable solutions represented by the solid 
lines and the unstable solutions represented by the dashed lines in both cases. The stable 
steady state solution, that the vehicle will acquire, is dependent on the initial conditions 
that are given to the vehicle. Using the results displayed in Figure 4-1 as an indication that 
xcg may also be a critical parameter, the subsequent task is to determine the relationship 


between the critical Froude number and xGz or f( Fn, x,,). 


B. BIFURCATION GRAPHS 
The first order of business іп the attempt to identify (Ри, хо») 1$ to simplify the 
exact pitchfork equation. This can be accomplished by replacing the cos@ and sin® terms 


in equation 4.3 with the following Taylor series expansions 


w lw 2ww -w 
іпб----------- 4.5 
IAN ди? >. 


и OU и. 
cos 21-—— 2 -————, 


- 4.6 
2 u 2H co 


which will allow the exact solution to be represented as follows; 


Z,( z,W —z4B)w! - 2u'C, A( M, — x ,Z; )w|w| 
-(2w ( Z, M, - M,Z,) * 2! Z,( 4$W -z,B))w -2w ((W - B)M (4.7) 
+(x,W — x,B)Z, - ((W - B)M; + (xW -x,B)Z,)w' 20 
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v 4 





1 
where C, now equals the terms = present in equation 4.3. Defining the following 


parameters 
аа 
Хов 7 Xo — Хв, 


ay 


^GB — 76 ^ 7p» 


г Й —z,B 2 z4W — 2,0, = Ô., 


х= —x,B-x.W-x0.250., 


ЖАРДЫ 


(4.8) 
(4.9) 
(4.10) 
(4.11) 
(4.12) 


(4.13) 


and substituting them into equation 4 7 will allow the coefficients in front of the various 


powers of w to be written as follows; 


w^ A 7 Z(z, M 2,95) 2 ZW - габ»), 


w|w|: 4, 2 2u'C5 A( M, - x ,Z, ), 


w: A, 2 -2 (i (Z, M; - M,Z,) t Z,(kwW — z,6,)), 
w^: A, 2 ди ( -6, M, - Z,( x,W — x40, )) 


w": А = (—Š,M, - Z,(x W — x,6,)). 
The resulting form of equation 4.7 can be written 


Ам? + Aw|w|+ Aw + À + Aw = 0, 


(4.14) 
(4.15) 
(4.16) 
(4.17) 
(4.18) 


(4.19) 


which is a generic pitchfork equation. If it is assumed that the submarine is neutrally 


buoyant, and equation 4.19 is divided by Z,z,,W, then the above coefficients can be 


redefined as; 


A, = 1, 


А, оС сое а вы 
ZW ga. 
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(4.20) 


(4.21) 





Á, 2 -2u/(—— +1), 
св ZW 
E Hu 2g. 
82св 
А, = са Ха 
8208 И 


Letting the critical parameter À be defined as 


^ 
- 


М 


3 
S765 





д = 


will allow equation 4.19 to be put into the following pitchfork solution form; 


и’ + YÀwlw| - 2i (1- QA )w - 2Bu^X  BAw^ 0, 
where, 


B-- Xgg& 


М 


c oe Ea 


° ZW | 


y 22uC,A( M; - x Z; е 


(4.22) 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


(4.27) 


(4.28) 


(4.29) 


To determine if the above manipulation and definitions are correct and valid, take the 


symmetric case that was developed in the previous chapter, i.e, x;,,-—0. This causes 


equation 4.26 to reduce to 


w(w^ 4 yAw|-2w(140X)) 20, 


(4.30) 


resulting in a solution always occurring at w=0, with two more solutions appearing at the 


bifurcation value 


1+CA, =0. 


Rearranging and substituting the appropriate values into equation 4.31 will yield 
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(4.31) 


” 


u Z; 
“GB M. Z, à а М, | 





(4.32) 


which is identical to equations 3.9 and 3.20, thus verifying the manipulation. 

To find the relationship between Fn and xcg the function A( w,À, B, and its partial 
derivative A (w, A, [), must be determined. These two functions must then be set equal to 
zero and to each other in order to obtain a function independent of w. As seen, equation 


4.30 is much too complex to complete these requirements, but by neglecting the terms 


Aw|w'| and Aw- in equation 4.26 a simplified pitchfork of the form, 
w^ -2w (1 QA )w - 2B0u^À = 0, (4.33) 
or in general terms Ai( w, A, B) 2 0, can be obtained which can be manipulated to meet the 
above mentioned requirements, thereby determining the critical (A,B) curve. By taking 
the partial derivative, h,,, of equation 4.33, the following equation results; 
h, 23w^ -2w (14 QA). (4.34) 


This equation can be set equal to zero, yielding 
" 
w — (4.35) 


which can be substituted into equation 4.33 to obtain the value of w where the function 


and it's partial derivative is equal to zero; 


384. 


= —  — ... .36 
К 2(1 СА) ў, 


This value of w can now be substituted back into equation 4.35 to obtain the desired 
Jf (&. B), or 

279 X 8 (1 СА)”. (4.37) 

Attempting to solve this equation analytically for A as a function of D, to determine 


the shape of this curve is difficult. Therefore, an approximation to this equation is needed 


for small values of B. By rearranging equation 4.37 into the form 
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: 4.38 
D D (4.38) 


an equation that determines B as a function of À is now produced. If this function is 
represented in a Taylor series expansion and simplified, then equation 4.38 can be written 


as 
е з кВ arc ))^, (4.39) 


which can be considered as an approximation to the (A,B) curve, with equation 4.37 as 
the exact (A,B) curve. Through inspection it can be seen that the general shape of this 
curve is that of a cusp. 

Figure 4-2 shows a comparison of four bifurcation curves. The curves represent the 
relationship between the critical speed U, which can be converted into the non dimensional 
parameter A using equation 4.25, and the value of xg, which also can be converted into a 
non dimensional parameter 8 using equation 4.27. The curves labeled exact bifurcation 
set, pitchfork bifurcation set, exact cusp and approximate cusp come from the numerical 
solution to equations 4.7, 4.26, 4.37 and 4.39 respectively. Observe that in the general 
vicinity of xGg=0 and near the critical speed the shape of each of these curves is indeed 
that of a cusp, with the peak occurring at the critical speed of 1.9 fps. Now, with the 
comparison conducted in Figure 4-2 complete, equation 4.7, the exact pitchfork solution 
will be used for the remainder of this analysis. 

Returning to equation 4.7, it can be seen that there is only one additional parameter, 
Cp, that the cusp curve produced using equation 4.39 does not incorporate (because zog 
апа хсв аге incorporated into the non dimensional parameters À and D). It is therefore 
necessary to conduct a sensitivity analysis to determine the effect that Cp has on the shape 
of the cusp curve. Figure 4-3 displays the results obtained when this sensitivity analysis 


was conducted. 
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Figure 4-3. Comparison of Exact Bifurcation Set for Different Values of Drag 
Coefficient. 
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As seen in the plot, the drag coefficient has no effect on the location of the peak of 
the cusp. The drag coefficient, however, does modify the slope of the cusp as displayed by 


the four cases shown in Figure 4-3. 


C. SOLUTION SETS 

[n this section, the equilibrium solutions to equations 4.3 and 4.5 will be investigated. 
As determined in the previous section, the critical parameters that control the location and 
shape of the pitchfork solution sets are A, В апа the drag coefficient C,. By running the 
FORTRAN program of Appendix A for several values of x,, Г the shape and trend of 
the solution sets can be determined. The result of this comparison, for a C, =0.0, is 
displayed in Figure 4-4. This plot supports the symmetry about x,, =0 displayed in the 
earlier cusp curves Observe that the critical or bifurcation point moves to a lower value of 
À as Хор increases, and that for low values of A the equilibrium solutions converge to the 
same value independent of the value of xg. 

When the comparison is made regarding the effect that the drag coefficient has on the 
solution sets, for a given value of xcg, Figure 4-5 is produced. This variation in drag 
coefficient for a constant value of xcg has a similar effect on the movement of the 
bifurcation point as fixing the drag coefficient and varying xcp. The difference is, however, 
evident by studying the stable solutions. The stable solutions for non-zero drag 
coefficients become more linear as the value of C, increases. This results from the 
increased effect of the w|w| term in equation 4.3 as the value of C, increases. 

Although Figure 4-5 shows the actual solution to the exact pitchfork equation it 1s 
not realistic since the effect of dive plane saturation is not considered. Once the dive 


planes saturate, the equations used to obtain the steady state solutions displayed in 
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Figure 4-4. Comparison of Exact Solution Set for Different Values of Xgb/L. 
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Figure 4-5. Comparison of Exact Solution Set for Different Values of Cp. 


62 


Figures 4-4 and 4-5 are no longer valid. Instead, equations 4.1 and 4.2 are modified, still 
remembering that the assumption of neutral buoyancy remains in affect, resulting in the 
following two equations; 


Z uw —C,Aw|w|+w°Z,8,,, = 0, (4.40) 


M uw — Cyx ,Aw|w| -Wx,, cos 0 — z,,W sin6 + M,u°d,,, = 0. (4.41) 
To obtain an equation that can be used to compute the steady state Ө solution while 
taking into consideration dive plane saturation requires solving equation 4.40 for w and 
then substituting this value into equation 4.41 which should yield an equation independent 
of w. The ability to perform this manipulation is hindered by the inclusion of the absolute 
value term in both of the above equations. However, if the absolute value terms are 
neglected by considering the case of Cp=0.0, then a single equation for the steady state 
value of O may be obtained. 
Neglecting the drag affects, and performing the algebra will yield an equation of the 
following form; 

( M,Z, - M,;Z, )w 8, * Z,W( x, cos0 * z,, sin8) = 0, (4.42) 
which may be used to determine the values of 0 once the dive planes have saturated. The 
validity of neglecting the drag terms is questionable, and an analysis that incorporates C 
will be conducted later in this section once the effect of varying xp is determined. Figure 
4-6 displays the changes that occur to the solution set of Figure 4-1 when saturation 15 
considered. Notice that saturation occurs when the solution set is at approximately nine 
degrees, and that when multiple solutions occur, the stable equilibrium states are 


represented by the saturation portions of the solution set. 
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Figure 4-6. Exact Solution Set for Xgb/L=-0.0001 With Saturation Included. 
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When the saturation condition is applied to the previously displayed solution sets for 
different values of xcp (Figure 4-4), the results shown in Figure 4-7 are produced. Items 
of interest in Figure 4-7 are that: 

1) The bifurcation point on each solution set moves to a lower value of A as 
the absolute value of xcg increases. 

2)  Therange where multiple solutions exist decreases as the absolute value of 
Xcp increases. 

3) | The amount of trim produced, in still water, by the given value of xcg can be 
obtained from the point where the upper an lower solution branches 
converge. 

Incorporating the saturation condition into the solution sets raises some additional 
concerns. To begin, what if the point of dive plane saturation occurs after the bifurcation 
points identified in Figure 4-4? If this condition did occur, referring to Figure 4-1 which 
displays the stable and unstable solution branches, there would be at least one value of À 
where there were [ms stable solution branches. To determine if this condition is possible 
a comparison of solution sets with and without saturation is necessary. This comparison is 
shown in Figure 4-8. As this plot depicts, the range between the bifurcation point and 
saturation point increases as the value оҒ хср increases. This result will allow the fore 
mentioned concern, of three stable solution branches, to be disregarded. 

A second concern is the validity of the previously developed bifurcation cusp curves 
Since the exact bifurcation equation, equation 4-7, does not include saturation, the cusp 
curves produced using this equation will not accurately predict the occurrence of multiple 
solutions. However, if equation 4.42 is expanded in a Taylor series, it can be manipulated 


into a form that will allow a cusp curve, which includes saturation, to be developed. 
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Figure 4-8. Comparison of Bifurcation Points and Dive Plane Saturation Points. 
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Performing this manipulation will produce the following equation; 


и" Б жн - —. (4.43) 
(Z;M, - MZ, )(9,,(1*(2,; 26)“ ) 

By converting 1/7 and xgg into the non dimensional parameters À and B the cusp 
curve shown in Figure 4-9 can be produced. This cusp displays the saturation point for 
both solution branches, and as seen it differs quite extensively from the exact bifurcation 
cusp without saturation. If Figure 4-9 is replotted to allow the peak area of the cusp to be 
more accurately represented, then Figure 4-10 is produced. By using these two cusp 


curves then the general form of the resulting solution set for any path through the cusp can 


be predicted It is this prediction that the subsequent section will address. 


D. PATH FORMULATION 

The ability to accurately predict how the vessel will respond to various changes in 
operating conditions is critical for proper control of the vessel. It is this prediction that 
necessitates the need for path analysis to be conducted. If the ballast condition on the 
vessel is held constant, xc;g --0.0001, and the speed of the vessel is varied, path A of 
Figure 4-11 is obtained. Neglecting for the moment the saturation cusp shown in 
Figure 4-11, it would be expected that a single steady state solution would exist for all 
values, of the parameter sgrt(X.), greater than the point where the path intersects the solid 
cusp curve. At the speed where the path intersects the cusp the bifurcation point occurs, 
and for speed values which place the path inside the cusp multiple steady state solutions 
occur. The results of this path analysis are depicted in Figure 4-1, presented earlier in this 
chapter. 

If attention is returned to the path through the saturation cusp of Figure 4-11, then 
the same type of predictions made for the path through the solid cusp can be made for the 


this cusp, with some slight differences. As the speed decreases the path intersects the 
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Figure 4-9. Bifurcation Cusp With Dive Plane Saturation Included. 
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Figure 4-10. Expanded Bifurcation Cusp With Saturation Included. 
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Figure 4-11. Path Formulation for Xgb/L=-0.0001. 
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upper portion of the dotted cusp at an approximate value of 1.08. At this value saturation 
occurs on one of the solution branches but there remains only one steady state solution 
since the path has not yet entered the cusp. As the speed continues to decrease, the path 
intersects and enters the saturation cusp at an approximate value of 0.99 with saturation of 
the other solution branch and the existence of multiple steady state solutions occurring. 
These results are evident in Figure 4-6. 

Now, if the speed of the vessel is held constant, and the loading conditions on the hull 
are changed, then the paths shown in Figure 4-12 are obtained. Path B, at a Froude 
Number of 1.1, is outside the cusp, therefore only a single steady state solution for each 
value of xcg L should exist. While path C, at a Froude Number of 1.0, passes through the 
cusp at + 0.0001 values of xcp ZL resulting in single steady state solutions for values less 
than or greater than + 0.0001, while values between + 0.0001 result in multiple steady 
state solutions. Through numerical computation the solution set as a function of xcpL is 
obtained with the results shown in Figure 4-13. 

As seen in Figure 4-13, path B, represented by the dotted line, does in fact result in a 
single steady state solution throughout the range of xcp/L, however, for path C this is not 
the case. Path C, represented by the solid line, exhibits multiple solutions between the 
Хсв І, values of — 0.0001, as predicted. This is a hysteresis curve, with the unstable 
solutions occurring between the corners of the curve. If the value of xog/L began 
at -0.0005 and progressed to the positive value of 0.0001, the solution set would remain 
on the curve with 0 taking on the values shown. Once the value of x5;g L becomes greater 
than 0.0001 the solution jumps from -9? to 7.5?. This same jump would be evident if the 
value of xcg L began positive and progressed negative. These jumps indicate that if the 


dive planes are held constant, then the vessel will instantaneously change its pitch angle. 
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Figure 4-12. Path Formulation at a Constant Speed With Xgb/L Varying. 
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Figure 4-13. Solution Set for Path Formulation. 
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To counter this effect the dive plane angle must be reversed. It is this dive plane reversal 
that the following chapter will discuss. 

Before proceeding to the next chapter it is necessary to determine the effect that drag 
has on the previously developed bifurcation cusps. Using equations 4.40 and 4.41, and 
substituting equation 4.4 into both equations will result in the following two equations; 

Z u? tan 0— C,Au* tan Atan A+ Zu Ò, =0, (4.44) 

M, uw - C,yx Au^ tan Gan 9 -WX cp COSO- Zc W ѕіпе+ М,и?5 = 0. (4.45) 

By using the MATLAB function fzeros contained in the program of Appendix A, the 

value of @ that solves equation 4.44 can be determined as a function of speed. This 

relationship can then be used to solve equation 4.45 to obtain the relationship between xcp 
апа и. 

With this relationship obtained the cusp curves of Figure 4-14 can be plotted It can 
be seen that for a given value of xcGgL increasing the drag coefficient forces the 
bifurcation point to a lower speed. An interesting item to note is the diamond shape peak 
on these cusps. This is a result of saturation not occurring, for a small range of speeds, in 
the small region around xgp/L=0. Path formulation through these cusps can be developed 


as shown earlier in this section. 
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Figure 4-14. Bifurcation Cusp With Saturation Included for Different Values of Cp. 
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V. DIVE PLANE REVERSAL 


A. STERN PLANE REVERSAL 

The ability of a submarine to avoid detection and perform its mission is dependent on 
its ability to maintain ordered depth. Modern asymmetric submarines exhibit an inherent 
phenomenon known as dive plane reversal, which is difficult to deal with. To introduce 
this physical occurrence, refer to Figure 5-1 [Ref. 14] and consider the case of wanting to 
dive the vehicle to a deeper depth. At moderate to high speeds the stern planes would be 
deflected by an angle 6,, this would in turn cause a force and moment to be developed due 
to the angle of attack of the planes. The vessel will respond by pitching downward, and 
will assume some angle of attack to the fluid flow causing hull forces and moments to 
develop. The stern plane and hull forces bring the ship to some steady state pitch angle 
which is obtained when the restoring moment and pitching moments balance. The vessel 
now has the ability to /7y itself to the ordered depth with the normal hull force assisting in 
pushing the vessel to a deeper depth. 

At low speed we would expect the same to be true, however it is not. If the stern 
planes are set to an angle 6,, this would cause a force and moment to be developed due to 
the angle of attack of the planes, the vessel will pitch downward creating a pressure field 
around the vessel due to the angle of attack with the fluid resulting in hull forces and 
moments to develop. The vessel once again steadies out at some steady state pitch angle 
when the restoring moment and pitching moments balance. Since the speed of the vessel is 


not as great as in the previous case, the hull forces and moments are not as large, therefore 
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Figure 5-1. Stern plane Reversal . 
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instead of the vessel flying itself to the ordered depth with the hull force assisting, the 
stern plane force actual pushes the vessel upward causing it to rise instead of dive, with 
some non-zero pitch angle. [Ref. 14] This condition was evident in the simulation shown 
in Figure 3.3. 

With this description of the vehicle response complete, it is obvious that to dive the 
vessel at low speeds, the opposite stern plane angle must be ordered. This will result in the 
vessel being pushed to the ordered depth at some bow up attitude. This required shift in 
dive plane angle is what is known as STERN PLANE REVERSAL. 

To tie together the previously completed bifurcation analysis with the physical 
occurrence just discussed, Figure 5-2 is produced. Figure 5-2 presents the pitch angle and 
stern plane angle required for straight and level flight with the bow planes centered. As the 
vessel's speed decreases, the metacentric moment ( M,0) makes it increasingly difficult to 
maintain a positive pitch angle on the hull. The stern plane angle is negative trying to 
produce a pitch-up moment, but the downward force of the stern planes requires an even 
larger pitch angle [Ref. 14]. Between the speeds of 1.08 and 1.18 knots, the required stern 
plane angle is beyond the normal operating range and the vessel cannot be flown straight 
and level. Below the critical speed, the stern plane and vessel pitch angle must reverse sign 
to maintain neutral trim. 

In a final effort to show the physical significance of stern plane reversal, the force 
input to the stern planes is plotted. If the control gains are allowed to change as the speed 
of the vessel changes, Figure 5-3 is produced. This plot shows the non-dimensional force 
per unit depth error input to the stern plane as a function of speed It is basically the 


control system gain used to set the stern plane angle. As can be seen, as the critical speed 
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Figure 5-2. Variation of Pitch and Stern plane Angles as a Function of Speed. 
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Figure 5-3. Control System Force Input to the Stern planes. 
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is approached the value of input tends to infinity requiring the stern planes to operate 
outside normal values. As the speed passes the critical value the sign of this input changes 
requiring the stern planes to operate in the opposite direction. This reversal in sign 15 
evident of the stern plane reversal phenomenon shown and discussed earlier in this 


chapter. 


B. BOW PLANE REVERSAL 

The reversal effect that was shown to occur with the stern planes can also be shown 
to occur with the bow planes. If the value of ax is allowed to tend to oo, indicating no stern 
plane control and only bow plane control, then the same analysis conducted throughout 
this thesis for stern plane reversal can be performed for bow plane reversal. If the critical 
Froude number, (sqrt(A)), is plotted as a function of a, then the results shown in 
Figure 5-4 are produced. This plot shows that for bow plane control only, the critical 
value approaches 2.05 as & approaches ©. 

To show that the same type of results can be obtained for the bow plane analysis, the 
case of neutral buoyancy, x¢g = 0 and dive plane saturation, with bow plane control only 
was investigated. Referring to Figure 5-5, it can be seen that at a Froude number of 2.05 
multiple steady state solutions occur. These results are similar to the results displayed in 
Figure 3-5. The major differences are that the critical point occurs at a Froude number of 
2.05 vice 1.05, and the maximum value of 0 is only + 4.5° vice + 9.5°. The reason for 


these changes is due to the values of Z, and M; being only one-half the respective value 


for the stern planes. 
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Figure 5-5. Steady State Values of 6 for Bow plane Control. 
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C. BIAS EFFECTS 

As a final area of analysis, the condition of operating near the surface, such as during 
periscope operations, must be investigated to determine the effects this will have on the 
previously developed bifurcation cusps. Operating near the surface will require the 
incorporation of a suction force and moment into the equations of motion. This force and 
moment is a function of the distance away from the surface, and the angle the vessel 
makes with the surface. Therefore, as a simple model of suction effects consider the case 
of operating with excess buoyancy. To conduct the analysis requires the use of equations 
4.1, 4.2 and 4.4. Performing the same substitution and solution technique as was used in 
Chapter V to obtain the saturation curves that incorporated the drag coefficient will result 
in solving the following equations to obtain the bifurcation cusps that would quasi-model 


surface effects; 


Zw - 7 pC, Au tan Oan 04 (W — B)cos0+ Za Ó,, = 0. (5.1) 


M. — - pCyyx, Au tan Өтап | - (x,W — x,B)cos0 (5.2) 
-(z,W — z,B)sin0+ Ma ó, , = 0 

If the value of Cp is set equal to zero then the saturation cusp produced in 

Figure 4-10 can be compared to the results obtained for several cases of excess buoyancy. 

Since it is the saturation cusp that controls the location of the bifurcation point, as shown 

in previous the chapter, the non saturation cusp will be neglected for the remainder of the 

analysis. Figure 5-6 shows the results obtained as excess buoyancy is increased. As can be 


seen, the cusp peak moves to down and to the right occurring at a lower speed and at a 


non zero value of xgg L. The body of the cusp also shifts to the right, as seen, with the 


85 


sqrt( Lambda) 


Ls 


о 
NO 
T mmm 


0.8 


BIASED BIFURCATION CUSP 


- for excess Биоуапсу 0.01 % of weight 


But | -- for excess buoyancy 0.02 % of weight 


-. for excess buoyancy 0.05 % of weight 





Xgb/L 
Figure 5-6. Bias Effects Caused by Operating Near the Surface. 
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x1 


entire cusp occurring in the positive region of the xgẹ L range for a value of 0.05 % 
excess buoyancy. 

To determine how the value of drag coefficient effects the cusp curves, consider the 
case of Cp = 0.3 shown in Figure 5-7. These cusps were produced for the same values of 
excess buoyancy used in Figure 5-6. Notice that the cusp curves shift to the right by the 
same value seen in Figure 5-6, but the peak value now occurs at an even a lower speed. 

Іп order to observe the sensitivity of the biased cusp to drag only, consider the case 
of setting the excess buoyancy to 0.01 % and varying the drag coefficient shown 
Figure 5-8. These results indicate that increased drag causes the critical speed to occur at 
a lower value. They also show that drag tends to counter the biasing effect that excess 
buoyancy was previously show to have. This is evident by the peak point drifting to the 
left as the value of drag coefficient increases. 

In summary, it can be stated that operating near the surface will tend to reduce the 
critical speed determined in previous chapters and the vessel will have to be placed in a 
different trim condition to counter surface effects, and that increasing drag also decreases 


the critical speed but tends to counter the trim condition required. 
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Figure 5-7. Bias Effects Caused by Operating Near the Surface Considering Drag. 
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Figure 5-8. Bias Effects Caused by Drag While Operating Near the Surface. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

e A method for analyzing submarine depth keeping and dive plane reversal has been 
developed and presented. This analysis identified the three phenomena, a positive real 
eigenvalue, a steady state pitchfork bifurcation, and an uncontrollable system, that lead to 
vessel instability. 

° The critical parameters have been identified that control the bifurcation. These 
parameters combine speed, loading conditions and hydrodynamic coefficients into 
non-dimensional values that may be applied to other models or vessels. 

е Since the non-dimensional Froude numbers are based on the metacentric height 
and LCG-LCB separation, Froude scaling can be used to apply this analysis to full scale 


vessels, allowing for the reduction of model testing during design phases. 


B. RECOMMENDATIONS 


Some recommendations for further research are as follows: 


e Develop and incorporate into the bifurcation analysis a more complete expression 
for the free-surface effects. 
e Conduct a Hopf bifurcation analysis to identify the critical parameters that lead to 


vehicle instability at high speeds. 
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APPENDIX A 


COMPUTER PROGRAMS 


% THIS PROGRAM IS THESIS1.M. IT SIMULATES THE RESPONSE OF THE 
DTRC SUBOFF 

% MODEL. IT PLACES THE POLES, CALCULATES THE CRITICAL SPEED, AND 
USES 

% PROGRAM THESISDE.M TO SIMULATE THE VERTICAL PLANE REPONSE. 


clg; 
global zg m zgb u U Zq Zw Zdlt Mq Mw B1 Mdlt K1 M bm xL Zval Mval rho 


yo ESTABLISH GEOMETRIC PARAMETERS 
i 2055602365.B1-—1556.2363; 

есі 92: 

р--32 2: 

m-W/g; 

rho=1.94: 

L=13.9792: 

xg=0;zg=0.1; 

zgb-0.1; 


% BEAM AND LENGTH DEFINITION FOR TRAPAZIODAL RULE USE 
bm=[0 485 658 778 871 945 1.01 1.06 1.18 1.41 1.57 1.66 1.67 1 67 1.67 1.63... 
1.37 919 448 195 188 168 132 053 0): 


x]2[0.1.2.3.4.5.6.71234 7.7143 10 15.1429 16 17 18 19 20 20.1... 
20.2 20.3 20.4 20.4167]; 

xl=(L/20)*x1; 

xL-xl-L/2; 


yo FACTORS TO CONVERT FROM NON-DIMENSIONAL VALUES TO 
DIMENSIONAL VALUES 

nd1-2.5*rho*L^2; 

пб гае сро” І. Э. 

nd3-2.5*rho*L^4, 

nd4-.5*rho*L^S; 


9] 


АІрһа=0; 


% NON-DIMENSIONAL HYDRODYNAMIC COEFFICIENTS 
Zqdnd--6.33e-4;Zwdnd--1.4529e-2;Zqnd-77.545e-3;Zwnd--1.391e-2; 
Zds-(-5.603e-3);Zdb-0.5*(-5.603e-3);Zdltnd-(Zds-*Alpha*Zdb); 
Madnd--8 8e-4;Mwdnd--5.61e-4;Mqnd--3.702e-3;Mwnd-1.0324e-2; 
Mds7(-0.002409),Mdb-0.5*(0.002409);Mdltnd-(Mds-Alpha*Mdb); 


% CONVERTION OF COEFFICENTS INTO DIMENSIONAL VALUES 
Zqd=nd3*Zqdnd;Zwd=nd2* Zwdnd;Zq=nd2*Zqnd;Zw=nd 1 *Zwnd, 
7ай-па1%7 айта; 


Maqd=nd4* Maqdnd:Mwd=nd3 *Mwdnd;Mq=nd3 *Mqnd;Mw=nd2*Mwnd; 
Mdlt=nd2*Mditnd; 


% FORMULATION OF A AND B MATRIX ELEMENTS 
Dv=(m-Zwd)*(Iy-Mqd)-(m*xg+Zqd)*(m*xg+Mwd); 

al 1 Dv=(ly-Mqd)*Zw+(m*xg+Zqd)*Mw; 

al 2Dv=(ly-Mqd)*(m+Zq)+(m*xg+Zqd)*(Mq-m* xg), 
al3Dv=-(m*xg+Zqd)* W; 
b1IDv=(Iy-Mad)*Zdlt+(m*xg+Zqd)*Madlt; 

a2 1Dv=(m-Zwd)* Mw+(m*xg+Mwd)* Zw; 
a22Dv-(m-Zwd)*(Mq-m*xg)*(m*xg*Mwd)*(m-*Zq); 
a23Dv--(m-Zwd)*W; 

b2Dv=(m-Zwd)* Mdlt+(m*xg+Mwd)* Zdlt; 


al l=al1Dv/Dv;al2=a12Dv/Dv;al3=a13Dv/Dv; 
а2 1=a2 1Dv/Dv,a22=a22Dv/Dv;a23=a23Dv/Dv; 
b1=b1 Dv/Dv:b2=b2Dv/Dv; 


% ESTABLISHMENT OF POLE LOCATIONS 
Еа 


% CALCULATION OF CRITICAL SPEED 
Ucr-sqrt(-(a13*b2-a23*bl1)*zgb/(all*b2-a21*bl)) 


% INPUT OF SPEED WHICH THE VESSEL WILL USE IN THE SIMULATION 
U-input('enter speed of vehicle ') 


% CALCULATION OF FROUDE NUMBER 
Fn-sqrt(U^2/(zgb*g)) 


25 


% CALCULATION OF MATRICIES FOR POLE PLACEMENT AND 
CONTROLABILLITY ANALYSIS 

и=5; 

a-[0 O 1 0;al3*zgb al1*u al12*u 0;a23*zgb a21*u a22*u O;... 

-u 1 0 0]: 

а1=[0 0 1 0;a13*zgb all *Uer al2*Ucr 0;a23*zgb a21*Ucr a22*Ucr 0... 
-Ucr 1 00]; 


b=[0;b1 *u%2;b2*u%2;0]; 
b1=[0;b1*Ucr^2;b2*Ucr^2;0],; 


КІ-ріасе(а,Б,р1); 

C=rank(ctrb(al,b1)); 

H=(Mw*Zdlt-Zw*Madlt)/(zgb*B 1 *Zdlt), 

theta l2acos(H*U^2)* 180/pi; 
Dss--(Zw/Zdlt)*tan(acos(H*U^2))*180/pi; 
Z-(Dss*Kl1(1)*acos(H*U^2)*K1(2)*U*tan(acos(H*U^2))).... 
(-K1(4)); 


% ESTABISHES VEHICLE MASS MATRIX 
M=[1 0 0 0;0 (m-Zwd) -Zad 0 ;0 -Mwd (Iy-Mad) 0 ;0 0 0 1): 


% SET TIME LIMITS AND INITIAL CONDITIONS FOR ODE SOLUTION 
t0=0;tf1=2500*L/U; 
x0=(0 0 0 1]; 


% SOLVES VEHICLE SIMULATION ODE'S 
[t1, X1]-ode4S('thesisde' tO,tf1,x0); 


% REDEFINE THE OUTPUT VECTOR ELEMENTS AND FORM DIVE PLANE 
RESPONSE 

м\м1=Х1(:,2);41=Х1(:,3);21=Х 1(:,4)461=Х1(:,1); 
d=-(K1(1)*thi+K1(2)*w1+K1(3)*q1+K1(4)*z1); 

11=length(d); 


for n=1:11, 

if d(n)> 4, 
d(n)=.4; 

elseif d(n)<-.4, 
а(п)--.4; 

end, 

end; 
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% PLOT THE RESULTS 

subplot(211); 

plot(t! *U/L,th] * 180/pi,t 1 *U/L,d* 180/pi,'--'):title"THETA/DELTA vs. TIME '): 
xlabel(‘non-dimensional time ');ylabel(‘theta/delta (degrees)'); grid; 
gtext('--delta');gtext({'for a speed of ' num2str(U) ' fps']); 


plot(tl *U/L.zl/L):title DEPTH vs. TIME);xlabel('non dimensional time"); 
vlabel('non dimensional depth ');grid; 

pause; 

subplot(111); 


% THIS IS PROGRAM THESISDE.M. IT CONTAINS THE SUBOFF MODEL 
SYSTEM OF 
% EQUATIONS TO BE SOLVED BY PROGRAM THESIS1.M 


% SET SYSTEM OF EQUATIONS FUNCTION 
function X]1dot=thesisde(tl,X1); 


% SET CONTROL LAW 
d1=-K1(1)*X1(1)-K1(2)*X1(2)-K1(3)*X1(3)-K1(4)*X 1(4); 


% APPLY DIVE PLANE SATURATION 
if d12.4, 
dl= 4; 
elseif dl <-.4, 
dl=- 4, 
end; 


% INCLUDE DRAG TERMS AND CALCULATE AREA AND CENTROID 
CD=0; 


if X1(2)==0, 
Х1(2)=1е-5, 

епа; 

if X1(3)==0, 
X1(3)=le-5; 

end, 


Zval-bm *((X1(2)-xL.*X1(3)).^3) /(abs(X1(2)-xL.*X1(3))); 
Mval-bm.*((( X1(2)-xL.*X1(3)).^3). *xL)./(abs(X1(2)-xL.*X1(3))); 


ZDragval-0;MDragval-0; 
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% trapaziodal integration 
for n=1:length(xL)-1, 


Zdragval=0.5*(Zval(n)+Zval(n+1))*(xL(n+1)-xL(n)); 
Mdragval=0.5 *(Mval(n)+Mval(n+1))*(xL(n+1)-xL(n)); 
ZDragval=ZDragval+Zdragval: 
MDragval=MDragval+Mdragval; 


end; 


Zdr=(0.5)*rho* CD*ZDragval; 
Mdr-(0.5)*rho*CD*MDragval; 


% VEHICLE SYSTEM OF EQUATIONS 

Fl2[X1(3); m*(U*X1(3)*zg*X1(3)^ 2)? Zq*U* X1(3)*Zw*U*X1(2)*U^2*Zdlt*d1-Zdr, 
m*(-zg*X1(2)*X1(3))*Mq*U*X1(3)* Mw*U*X1(2)- 

Bl*zgb*sin(X1(1))*U^2*Mdlt*dl -Mdr; 
-U*sin(X1(1))*X1(2)*cos(X1(1))]; 


mpr-inv(M ); 
XIdot-mpr*F 1; 
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PROGRAM STEADY 


STEADY STATE SOLUTIONS OF SUBMARINES 
IN THE DIVE PLANE AT LOW SPEEDS 


INPUT DATA 


ICON = 1 : COMPUTATION OF SOLUTION SETS (S.S) 
2: COMPUTATION OF BIFURCATION GRAPHS (B.G ) 


IVAR = 1: Fn VARIATION 
2: Xgb VARIATION 


REAL MW,LENGTH,MASS,MDS,MDB,MD,LAMBDA 
DIMENSION B(25),X(25),BX(25) 
GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 


RHO =1.94 

LENGTH= 13.9792 

WRITE (*,*)' ENTER CDZ' 

READ (*,*) CDZ 

CDZ - CDZ*0.5*RHO 

BUO- 1556.2363 

WRITE(*,*) ENTER VALUE OF EXCESS BUO' 
READ (*,*) ЕАС 

WEIGHT-FAC*BUO 

MASS = WEIGHT/32.2 

XB =00 

ZB =0.0 

WRITE (*,*) 'ENTER THE METACENTRIC HEIGHT ZGB' 
READ (*,*) ZG 

ZW =-0.013910*0.5*RHO*LENGTH**2 
ZDS --0.005603*0.5*RHO*LENGTH**2 
ZDB --0.005603*025*RHO*LENGTH**2 
MW =0.010324*0.5*RHO*LENGTH**3 
MDS =-0.002409*0.5*RHO*LENGTH**3 
MDB = 0.002409*0 25*RHO*LENGTH**3 


WRITE (*,*) ENTER THE VALUE OF ALPHA 


READ (*,*) ALPHA 
ZD=ZDS+ALPHA*ZDB 
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MD=MDS+ALPHA*MDB 
С 
С DEFINE THE LENGTH X AND BREADTH B TERMS FOR THE 
INTEGRATION 
с 
Х( 1)-20.4167 
Х( 2)-20.4000 
Х( 3)=20.3 
Х( 4)-20.2 
Х( 5)-20.1 
X( 6)-20.0 
Х(7)-19.0 
Х( 8)-18.0 
Х( 9)-17.0 
X(10)=16.0 
Х(11)=15.1429 
Х(12)=10.0 
Х(13)-7.7143 
Х(14)-4.0 
Х(15)-3.0 
Х(16)-2.0 
Х(17)-1.0 
Х(18)-0.7 
Х(19)-0.6 
X(20)=0.5 
X(21)=0.4 
X(22)=0.3 
X(23)=0.2 
X(24)=0.1 
X(25)=0.0 


B( 1)=0.00000 
B( 2)=0.03178 
В( 3)-0.07920 
В( 4)-0.10074 
В( 5)-0.11243 
В( 6)-0.11724 
В( 7)-0.26835 
В( 8)-0.55025 
В( 9)-0.81910 
В(10)=0.97598 
В(11)-1.00000 
В(12)-1.00000 
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В(13)=1.00000 
В(14)-0.99282 
B(15)=0.94066 
B(16)=0.84713 
B(17)=0.70744 
B(18)=0.63514 
B(19)=0.60352 
B(20)=0.56627 
B(21)=0.52147 
B(22)=0.46600 
B(23)=0.39396 
B(24)=0 29058 
B(25)=0.00000 


DO 1 I=1,25 
B(1)=B(I)* 1.6667 
X(1)=(-X(1)+10.0)*_LENGTH/20.0 
ВХ(1)-В(1)“Х(1) 

1 CONTINUE 

CALL TRAP(18,B,X,AREA) 
CALL TRAP(18,BX,X,CH) 
XA=CH/AREA 
WRITE(*,*) 'AREA EQUALS | AREA 
WRITE(*,*) XA EQUALS „ХА 
IFC=0 


OPEN RESULTS FILES 
OPEN (10,FILE='R10.RES'STATUS="NEW’) 
OPEN (11,FILE='R11.RES'SSTATUS='NEW’) 
OPEN (12,FILE='R12.RES',STATUS='NEW’) 
OPEN (13,FILE='R13.RES'SSTATUS="NEW’) 
OPEN (14,FILE='R14.RES' STATUS="NEW’) 
OPEN (15,FILE='R15 RES STATUS="NEW)) 
OPEN (16,FILE='R16.RES'STATUS="NEW’) 
OPEN (20,FILE='C20.RES' STATUS="NEW’) 
OPEN (21,FILE='S21.RES' STATUS="NEW’) 
OPEN (22,FILE='S22. RES'STATUS="NEW’) 
OPEN (23,FILE='S23 RES'STATUS="NEW)) 
OPEN (24,FILE='S24.RES'STATUS="NEW’) 
OPEN (31,FILE-'S31.RES', STATUS-'NEW)) 
OPEN (32,FILE='S32.RES' STATUS="NEW’) 
OPEN (33,FILE='S33.RES'STATUS="NEW’) 
OPEN (34,FILE='S34.RES' STATUS="NEW’) 
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СИРО) 


OPEN (41,FILE='C41.RES' STATUS=NEW’) 
OPEN (42, FILE='C42.RES' STATUS='NEW’) 
OPEN (43, FILE='C43.RES' STATUS='NEW’) 


READ DATA INTERACTIVELY 


WRITE (*,1001) 

READ (*,*) ICON 

IF (ICON.EQ 2) GO TO 331 
WRITE (*,1009) 

READ (*,*) IVAR 

GO TO (331,332), IVAR 


331 WRITE (*,1002) 


READ (*,*) ОМІ 
WRITE (*,1003) 
READ (*,*) UMAX 
WRITE (*, 1004) 
READ (*,*) ITER 
JU=ITER 

IF (ICON EQ 2) GO TO 332 
WRITE (*,1010) 
READ (*,*) XG 
XG=XG*LENGTH 
GO TO 300 


332 WRITE (*,1005) 


READ (*,*) XGMIN 
WRITE (*,1006) 

READ (*,*) XGMAX 
WRITE (*,1004) 

READ (*,*) ITER 
JXG-ITER 
XGMIN-XGMIN*LENGTH 
XGMAX-XGMAX*LENGTH 
IF (ICON.EQ.2) GO TO 334 
WRITE (*,1012) 

READ (**) U 


334 CONTINUE 


GO TO 300 


300 WRITE (*,1013) 


READ (*,*) ISET 
IF (ICON.EQ 2) GO TO 370 
GO TO (350,360,400), ISET 
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С 
C EXACT SOLUTION SET 
C 
350 DO 2 I=1,ITER 
WRITE (*,2000) LITER 
GO TO (351,352), IVAR 
U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 
GO TO 355 
XG=KGMIN+(XGMAX-XGMIN)*(I-1)/(ITER-1) 
355 Al= ZW*MD-MW*ZD 
A2=-CDZ* AREA*(MD-XA*ZD) 
A3=-(XG-XB)*WEIGHT*ZD 
A4= (ZG-ZB)*WEIGHT*ZD 
CALL 
EXSOLS(IVAR.LL,IFC,U,XG,A1,A2,A3,A4,CDZ, AREA,ZW,ZD,ZG,MW, 
& MD.BUO,XA) 
IF (LL.GT.1) IFC=1 
2 CONTINUE 
СО ТО 5000 


>> 
ол 
— 


+2 
Un л 
г 


С 
C PITCHFORK SOLUTION SET 
С 
360 DO 3 I=1,ITER 
WRITE (*,2000) LITER 
GO TO (361,362), IVAR 
361 U-UMIN-*(UMAX-UMIN)*(I-1Y(ITER-1) 
GO TO 365 
362 XG-XGMIN-*(XGMAX-XGMIN)*(I-1)(ITER-1) 
365 Al- ZD*(ZG-ZB)*WEIGHT 
А2= 2 0*U*U*U*CDZ* AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2 0*U*U*ZD*(XG-XB)*WEIGHT 
AS--ZD*(XG-XB)* WEIGHT 
CALL PITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,A5) 
IF (LL.GT 1) IFC=1 
3 CONTINUE 
GO TO 5000 
е 
C SIMPLE PITCHFORK SOLUTION SET 
С 
400 DO 401 I=1,ITER 
WRITE (*,2000) LITER 
GO TO (402,403), IVAR 
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402 U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 


403 
405 


GO TO 405 


XG=XGMIN+(XGMAX-XGMIN)*(I-1)/(ITER-1) 
Al= ZD*(ZG-ZB)*WEIGHT 


A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2.0*U*U*ZD*(XG-XB)*WEIGHT 

CALL SIMSLS(IVAR,LL,IFC,U, XG,A1,A3,A4) 

IF (LL.GT.1) IFC21 


401 CONTINUE 
GO TO 5000 
370 GO TO (380,390,410,430), ISET 


C 


C EXACT BIFURCATION SET 


С 


380 DO A I-1.JXG 


WRITE (*,2000) LJXG 
XG=XGMIN+(XGMAX-XGMIN)*(I-1)/(JXG-1) 
ІН5-0 

DO 5 J=1,JU 


U=UMIN+(UMAX-UMIN)*(J-1)/(JU-1) 
Al= ZW*MD-MW*ZD 
A2=-CDZ* AREA*(MD-XA*ZD) 
A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4= (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM,A1,A2,A3,A4,U) 
IF (J. NE.1) GO TO 381 
UO-U 
NUMO=NUM 
GO TO 5 
UN-U 
NUMN-NUM 
NDIF-IABS(NUMO-NUMN) 
IF (NDIF.EQ 2) GO TO 382 
UO-UN 
NUMO-NUMN 
GO TO 5 
U1-UO 
U2-UN 
DO 61 JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al= ZW*MD-MW*ZD 
A2=-CDZ* AREA*(MD-XA*ZD) 
A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
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А4- (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM, A1,A2,A3,A4,U) 
IF (JJ NE.1) GO TO 383 
UO=U 
NUMO=NUM 
GO TO 61 
383 | UN-U 
NUMN=NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF.EQ 2) GO TO 384 
UO-UN 
NUMO-NUMN 
61 CONTINUE 
GO TO 5 
384  UI-UO 
U2-UN 
DO 62 JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/JU-1) 
Al= ZW*MD-MW*ZD 
A2=-CDZ* AREA*(MD-XA*ZD) 
A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4= (ZG*WEIGHT-ZB*BUO)*ZD 
CALL EXNUMS(NUM,A1,A2,A3,A4,U) 
IF (JJ.NE.1) GO TO 386 
UO-U 
NUMO=NUM 
GO TO 62 
386 UN=U 
NUMN=NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF EQ 2) GO TO 385 
UO=UN 
NUMO=NUMN 
62 CONTINUE 
GOTOS 
385 UP=(UO+UN)/2.0 
[HS=IHS+1 
UPL=SQRT(UP/(32.2*LENGTH)) 
XGL=XG/LENGTH 
IPF=10+IHS 
WRITE (IPF, 100) XGL,UP 
UO=UN 
NUMO=NUMN 
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5 CONTINUE 
4 CONTINUE 
GO TO 5000 


Ф 


C PITCHFORK BIFURCATION SET 


С 


390 DO6I-1,JXG 


391 


092 


WRITE (*,2000) LJXG 
XG-XGMIN-(XGMAX-XGMIN)*(I-1)/(JXG-1) 
IHS-0 
DO 7 J=1,JU 
U=UMIN+(UMAX-UMIN)*(J-1)/(JU-1) 
Al- ZD*(ZG-ZB)* WEIGHT 
A2- 2.0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3--2*U*U*(U*U*(ZW*MD-MW*ZD)*ZD*(ZG-ZB)* WEIGHT) 
A47 2.0*U*U*ZD*(XG-XB)* WEIGHT 
AS--ZD*(XG-XB)* WEIGHT 
CALL PINUMS(NUM,A1,A2,A3,A4,A5,U) 
IF (J.NE.1) GO TO 391 
UO-U 
NUMO=NUM 
GO TO 7 
UN-U 
NUMN-NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF.EQ.2) GO TO 392 
UO-UN 
NUMO=NUMN 
GO TO 7 
U1-UO 
U2-UN 
DO 71 JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al= ZD*(ZG-ZB)*WEIGHT 
A2= 2.0*U*U*U*CDZ* AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4= 2.0*U*U*ZD*(XG-XB)*WEIGHT 
AS--ZD*(XG-XB)*WEIGHT 
CALL PINUMS(NUM,A1,A2,A3,A4,AS,U) 
IF (JJ. NE.1) GO TO 393 
UO-U 
NUMO-NUM 
GO TO 71 
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393 | UN-U 
NUMN-NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF.EQ.2) GO TO 394 
UO=UN 
NUMO=NUMN 
71. CONTINUE 
STOP 1111 
394 Ul=UO 
U2=UN 
DO 72 JJ=1,JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 
Al- ZD*(ZG-ZB)*WEIGHT 
A2- 2.0*U*U*U*CDZ* AREA*(MD-XA*ZD) 
A3--2*U*U*(U*U*(ZW*MD-MW*ZD)*ZD*(ZG-ZB)* WEIGHT) 
A4= 2.0*U*U*ZD*(XG-XB)*WEIGHT 
AS--ZD*(XG-XB)* WEIGHT 
CALL PINUMS(NUM,A1,A2,A3,A4,AS,U) 
IF (JJ NE 1) GO TO 396 
UO=U 
NUMO=NUM 
GO TO 72 
396 | UN-U 
NUMN-NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF.EQ 2) GO TO 395 
UO-UN 
NUMO-NUMN 
72. CONTINUE 
STOP 1112 
395 UP=(UO+UN)/2.0 
IHS-IHS-*1 
UPL-SQRT(UP/(32.2*LENGTH)) 
XGL-XG/LENGTH 
IPF=10+IHS 
WRITE (IPF,100) XGL,UP 
UO=UN 
NUMO=NUMN 
7 CONTINUE 
6 CONTINUE 
GO TO 5000 
C 
C SIMPLE PITCHFORK BIFURCATION SET 
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С 


410 DO 411 I-1,JXG 


413 


414 


416 


415 


417 


WRITE (*,2000) LJXG 
XG=XGMIN+(XGMAX-XGMIN)*(I-1)/(JXG-1) 
ІН5-0 

DO 412 J=1,JU 


U=UMIN+(UMAX-UMIN)*(J-1)/(JU-1) 
Al= ZD*(ZG-ZB)* WEIGHT 
A3--2*U*U*(U*U*(ZW*MD-MW*ZD)*ZD*(ZG-ZB)* WEIGHT) 
A4- 2.0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SINUMS(NUM, A1,A3,A4,U) 
IF (J.NE.1) GO TO 413 
UO=U 
NUMO=NUM 
GO TO 412 
UN=U 
NUMN=NUM 
NDIF=IABS(NUMO-NUMN) 
IF (NDIF.EQ.2) GO TO 414 
UO=UN 
NUMO=NUMN 
GO TO 412 
U1-UO 
U2-UN - 
DO 415 JJ=1,JU 
U-UI-(U2-U1)*(JJ-1)/(JU-1) 
Al- ZD*(ZG-ZB)*WEIGHT 
A3--2*U*U*(U*U*(ZW*MD-MW*ZD)*ZD*(ZG-ZB)* WEIGHT) 
A4- 2.0*U*U*ZD*(XG-XB)* WEIGHT 
CALL SINUMS(NUM,A1,A3,A4,U) 
IF (JJ.NE.1) GO TO 416 
UO-U 
NUMO-NUM 
GO TO 415 
UN-U 
NUMN-NUM 
NDIF-IABS(NUMO-NUMN) 
IF (NDIF EQ 2) GO ТО 417 
UO-UN 
NUMO-NUMN 


: CONTINUE 


$ТОР 1111 
=) 
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U2-UN 
DO 418 JJ=1,JU 
U-UI-*(U2-U1)*(JJ-1)/(JU-1) 
Al- ZD*(ZG-ZB)* WEIGHT 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)* WEIGHT) 
A4= 2. 0*U*U*ZD*(XG-XB)*WEIGHT 
CALL SINUMS(NUM.AI,A3,A4,U) 
IF (JJ NE.1) GO TO 419 
UO=U 
NUMO=NUM 
GO TO 418 
419 | UN-U 
NUMN-NUM 
NDIF-IABS(NUMO-NUMN) 
IF (NDIF.EQ.2) GO TO 420 
UO-UN 
NUMO=NUMN 
418 CONTINUE 
STOP 1112 
420 UP=(UO+UN)/2.0 
IHS=IHS+1 
UPL=SQRT(UP/(32.2*LENGTH)) 
XGL=XG/LENGTH 


IPF=10+IHS 
WRITE (IPF,100) XGL,UP 
UO=UN 
NUMO=NUMN 
412 CONTINUE 
411 CONTINUE 
GO TO 5000 
C 
С ANALYTIC BIFURCATION SET 
C 


430 WRITE (*,1017) 
READ (*,*) ICUSP 
GO TO (432,433), ICUSP 
Ç 
C "EXACT" CUSP 
C 
432 DO 431 I-1.JU 
- WRITE (*,2000) LJU 
U=UMIN+(UMAX-UMIN)*(I-1)/(JU-1) 
DELTA=32.2*(ZW*MD-MW*ZD)/(ZD* WEIGHT) 
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LAMBDA-U*U/(2.2*(ZG-ZB)) 
BETA2=(8.0*U*U*(1.0+DELTA*LAMBDA)**3)/(27.0*LAMBDA **2) 
IF (BETA2.LT.0.0) GO TO 431 

ВЕТА-5ОЕТ(ВЕТА2) 

XGB-U*U*BETA/32.2 

WRITE (10,100) U,XGB/LENGTH 


43] CONTINUE 
GO TO 5000 
C 
С "APPROXIMATE" CUSP 
e 


433 DO 434 I=1,JU 
WRITE (*,2000) LJU 
U-UMIN-*(UMAX-UMIN)*(I-1)/JU-1) 
DELTA-32 2*(ZW*MD-MW*ZDY(ZD* WEIGHT) 
LAMBDA-U*U/(32 2*(ZG-ZB)) 
BET A2=(8.0*U*U*(1.0+DELTA*LAMBDA)**3)/27.0 
IF (BETA2.LT.0.0) GO TO 434 
ВЕТА-5ОКТ(ВЕТА2) 
XGB=U*U*BETA/32.2 
WRITE (10,100) U,XGB/LENGTH 
434 CONTINUE 
C 
5000 STOP 
1001 FORMAT (‘ENTER 1: SOLUTION SETS — '/ 
| '  2:BIFURCATION GRAPHS ’) 
1009 FORMAT (' ENTER 1 :U VARIATION '/ 
1 ' 2:XGVARIATION ) 
1002 FORMAT (' ENTER MINIMUM VALUE OF U’) 
1003 FORMAT ('‘ ENTER MAXIMUM VALUE OF U’) 
1004 FORMAT (' ENTER NUMBER OF INCREMENTS’) 
1005 FORMAT (' ENTER MINIMUM VALUE OF XG/L) 
1006 FORMAT (' ENTER MAXIMUM VALUE OF XG/L) 
1010 FORMAT (‘ENTER VALUE OF XGIL’) 
1012 FORMAT (' ENTER VALUE OF U’) 
1013 FORMAT (ENTER 1: EXACT COMPUTATION — '/ 
| '  2:PITCHFORK APPROXIMATION'/ 
2 '  3:SIMPLE PITCHFORK  '/ 
3 '  4:ANALYTIC SET ' 
1017 FORMAT (' ENTER 1 : EXACT CUSP '/ 
| '  2:APPROXIMATE CUSP ') 
2000 FORMAT (215) 
100 FORMAT (2E15.5) 
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ве) 


Ga СЭ 


SUBROUTINE TRAP(N,A,B,OUT) 


NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 


DIMENSION A(1),B(1) 
NI=N-1 
OUT=0 0 
DO I I-L,NI 
OUT 1=0.5*(A(1)+A(I+1))*(B(I+1)-B(D) 
OUT -OUT-OUTI 
| CONTINUE 
RETURN 
END 


SUBROUTINE EXSOLS(IVAR,L,IFC,U,XG,AI,A2,A3,A4,CDZ,AREA, ZW, 
&ZD.ZG,MW,MD,BUO, XA) 


IT COMPUTES AND PRINTS THE EXACT SOLUTION SET THETA 
(DEG.) VERSUS U OR XG 


REAL MW,MD,PRNT,PI 
DIMENSION VF(5,2) 

Р1=4 0*АТАМ(1.0) 

IF (IVAR EQ 1) PRNT-U 
IF (IVAR EQ 2) PRNT=XG 


FIND FIRST ESTIMATE OF SOLUTIONS 


L=0 
VA=-89 5 
VA-VA*PI/180.0 
VAO-VA 
VO-THETEQ(1,VA,A1,A2,A3, A4,U) 
VAMIN--89.5 
VAMAX=+89.5 
IV A=5000 
DO 10 I=2,IVA 
AI-I 
VA-VAMIN*(VAMAX-VAMIN)*(I-1)/(IVA-I) 
VA=AI-90.5 
VA-VA*PI/180.0 
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9 


УАМ-УА 
VN-THETEQ(1,VA,A1,A2,A3, A4,U) 
VP-VO*VN 

IF (VP.GE.0.0) GO TO 11 

L=L+1 

VF(L,1)=VAO 

VF(L,2)=VAN 


VO=VN 


VAO=VAN 


10 CONTINUE 


C EXACT COMPUTATION OF SOLUTIONS VIA NEWTON'S METHOD 


G 


С 


Е=1.Е-8 
IEND-500 
DO 20 J=1,L 


X=(VF(J,1)+VF(J,2))/2.0 


X1=X 
IXX=0 
IF (IXX EQ.0) GO TO 35 


F=THETEQ(1,X,A1,A2,A3,A4,U) 
FDER=THETEQ(2,X,A1,A2,A3,A4,U) 
DO 30 K=1,IEND 
IF (FDER.EQ.0.0) STOP 1001 
DX=F/FDER 
X1=X-DX 
F=THETEQ(1,X1,A1,A2,A3,A4,U) 
FDER=THETEQ(2,X1,A1,A2,A3,A4,U) 
IF (F.EQ.0.0) GO TO 35 
A-ABS(XI-X) 
IF (A-E) 35,35,40 
Х-ХІ 


30 CONTINUE 


GO TO 20 


35 SOL-XI*180/PI 
WRITE (*,*) X1 


JPR-104J 

WD=U*TAN(X1) 

DLT=((CDZ* AREA*WD* ABS(WD))-(ZW*U*WD))/(ZD*U*U) 
IF ((ABS(DLT).GT.0.4). AND.(CDZ.NE.0.0)) THEN 
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CALL 
SATUR(U,XG,CDZ, AREA,ZW,ZD,ZG,MW,MD,BUO, XA,DLT,PRNT) 
ELSEIF ((ABS(DLT).GT.0.4). AND.(CDZ.EQ.0.0)) THEN 
CALL 
SATURI(U,XG,CDZ, AREA,ZW,ZD,ZG,MW,MD,BUO, XA DLT PRNT) 
ENDIF 
IF (L.EQ.1. AND.IFC EQ.0) WRITE (10,100) PRNT,SOL,DLT*180/PI 
IF (L.EQ. 1. AND IFC.EQ.1) WRITE (16,100) PRNT,SOL,DLT*180/PI 
IF (L.GT.1) WRITE (JPR,100) PRNT,SOL,DLT*180/PI 
20 CONTINUE 
RETURN 
100 FORMAT (3E15.5) 
END 


SUBROUTINE 

SATUR(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 

REAL MW,MD,PRNT 
C 
C USED TO COMPUTE SATURATION CONDITIONS FOR THE DIVE PLANE 
WITH DRAG 
С COEFFICENT INCLUDED 
(С 

CD=CDZ 

IF (DLT GT.0.0) THEN 
W1=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD* AREA)*ZD*U*U*0.4)) 

& /(2*(-CD*AREA)) 

IF (W1.GT.0.0) THEN 

A=MW*U*W1-CD*XA* AREA*W1*W1-XG*BUO+MD*U*U*0.4 
=-2.0*ZG*BUO 
C=MW*U*W1-CD*XA* AREA*W1*W1+XG*BUO+MD*U*U*0.4 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL1=2.0* ATAN(X1) 
X2=(-B-SORT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0* ATAN(X2) 
WRITE (21,200) PRNT,SOL1,SOL2 
ENDIF 
W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD* AREA) *ZD*U*U*0.4)) 
& /2*(-CD* AREA)) 

IF (W2.GT.0.0) THEN 
A-MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO*MD*U*U*0.4 
B=-2.0*ZG*BUO 
C=MW*U*W2-CD*XA* AREA*W2*W2+XG*BUO+MD*U*U*0.4 
X1=(-B+SOQRT((B1*B1)-4*(A*C)))/(2.0*A) 
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SOL1=2.0*ATAN(X1) 
X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0* ATAN(X2) 
WRITE (22,200) PRNT,SOL1,SOL2 
ENDIF 
W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0.4)) 
& /(2*(CD*AREA)) 

IF (W3.LT.0.0) THEN 
A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*0.4 
B=-2.0*ZG*BUO 
C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*0.4 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL1=2.0*ATAN(X]1) 
X2-(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0* ATAN(X2) 
WRITE (23,200) PRNT,SOL1,SOL2 
ENDIF 
W4=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(CD* AREA)*ZD*U*U*0.4)) 
& /(2*(CD*AREA)) 
IF (W4.LT.0.0) THEN 
A=MW*U*W4+CD*XA* AREA*W4*W4-XG*BUO+MD*U*U*0.4 
--2 0*ZG*BUO 
C=MW*U*W4+CD*XA *AREA*W4*W4+XG*BUO+MD*U*U*0.4 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL1=2.0* ATAN(X1) 
X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0*ATAN(X2) 
WRITE (24,200) PRNT,SOLI,SOL2 
ENDIF 
ELSE 
WI-(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 
& (-0.4)))/(2*(-CD* AREA)) 

IF (W1.GT.0.0) THEN 
=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*(-0.4) 
--2.0*ZG*BUO 

C-MW*U*WI-CD*XA*AREA*WI*WI-XG*BUO*MD*U*U*(-0 4) 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL1=2.0* ATAN(X1) 

X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL2=2.0* ATAN(X2) 

WRITE (31,200) PRNT,SOL1,SOL2 

ENDIF 

W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD* AREA)*ZD*U*U* 
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& (-0.4))) (2*(-CD* AREA)) 

IF (W2.LT.0.0) THEN 
A=MW*U*W2-CD* XA *AREA*W2*W2-XG*BUO+MD*U*U*(-0.4) 
B=-2.0*ZG*BUO 
C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*(-0.4) 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 

SOL 1=2.0*ATAN(X1) 
X2-(-B-SORT((B1*B1)-4*(A*C))/(2.0* A) 
SOL2-2.0*ATAN(X2) 
WRITE (32,200) PRNT,SOL21,SOL22 
ENDIF 
W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD* AREA)*ZD*U*U* 
& (-0.4)))/(2*(CD* AREA)) 
IF (W3.LT.0.0) THEN 
A=MW*U*W3+CD*XA* AREA*W3*W3-XG*BUO+MD*U*U*(-0.4) 
=-2.0*ZG*BUO 
C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*(- 
04) 
X1=(-B+SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL1=2.0*ATAN(X]1) 
X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0* ATAN(X2) 
WRITE (33,200) PRNT,SOL1,SOL2 
ENDIF 
W4=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD* AREA)*ZD*U*U* 
& (-0.4)))/(2*(CD* AREA)) 

IF (W4 LT.0.0) THEN 
A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*(-0.4) 
B=-2.0*ZG*BUO 
C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*(- 

0.4) 
X1=(-B+SOQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL1=2.0*ATAN(X]1) 
X2=(-B-SQRT((B1*B1)-4*(A*C)))/(2.0*A) 
SOL2=2.0* ATAN(X2) 
WRITE (34,200) PRNT,SOL1,SOL2 

ENDIF 

ENDIF 
200 FORMAT (3E15.5) 
END 


SUBROUTINE 
SATURI(U,XG,CDZ, AREA.ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 
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OOC) 


Су С) СУС 


eC) 


REAL MW,MD,PRNT 
COMPUTES DIVE PLANE SATURATION NEGLECTING DRAG 


Р1=4.0*АТАМ(1.0) 

A1-MD*U*U*0.4-MW*U*U*ZD*0.4-XG*ZW*BUO 
=-MD*U*U*0.4+MW*U*U*ZD*0.4-XG*ZW*BUO 

B1--2.0*ZG*ZW*BUO 

C1--MW*U*U*ZD*0.4-XG*ZW*BUO*MD*U*U*ZW*0.4 

C2-MW*U*U*ZD*0.4-XG*ZW*BUO-MD*U*U*ZW*0.4 

XCD1=((-B1/A1)+SQRT(((B1/A1)**2)-4.0*(C1/A1)))/2.0 

XCD2=((-B1/A1)-SQRT(((B1/A1)**2)-4.0*(C1/A1)))/2.0 

XCD3=((-B1/A2)+SQRT(((B1/A2)**2)-4.0*(C2/A2)))/2.0 

XCD4=((-B1/A2)-SQRT(((B1/A2)**2)-4.0*(C2/A2)))/2.0 

XC1=2.0* ATAN(XCD1)*180/PI 

XC2-2.0* ATAN(XCD2)*180/PI 

XC3-2 0*ATAN(XCD3)*180/PI 

XC4-2.0* ATAN(XCD4)*180/PI 

WRITE (20,100) PRNT,XCI,DLT*180/PI 

WRITE (41,100) PRNT,XC2,DLT*180/PI 

WRITE (42,100) PRNT,XC3,DLT* 180/PI 

WRITE (43,100) PRNT,XC4,DLT*180/PI 


100 FORMAT (3E15.5) 


END 


SUBROUTINE PITSLS(IVAR;LL,IFC,U,XG,A1,A2,A3,A4,A5) 


IT COMPUTES AND PRINTS THE PITCHFORK SOLUTION SET THETA 
(DEG.) VERSUS U OR XG 


DIMENSION VF(5,2) 
PI-4.0* ATAN(1.0) 

IF (IVAR.EQ.1) PRNT-U 
IF (IVAR.EQ.2) PRNT-XG 


FIND FIRST ESTIMATE OF SOLUTIONS 


L-0 
VA--89.5 

VA-VA*PI/180.0 

VAO-VA 
VO=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 
DO 10 I=2,180 
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AI-I 
VA-AI-90.5 
VA=VA*PI/180.0 
VAN=VA 
VN-PITCEQ(1,VA,A1,A2,A3,A4,AS,U) 
VP-VO*VN 
IF (VP GE.0.0) GO TO 11 
L=L+1 
VF(L,1)=VAO 
VF(L,2)=VAN 
11 VO-VN 
VAO-VAN 
10 CONTINUE 
Є 
C EXACT COMPUTATION OF SOLUTIONS VIA NEWTON'S METHOD 
С 
Е=1 Е-5 
IEND-500 
DO 20 J=1,L 
X=(VF(J,1)+VF(J,2))/2.0 
F-PITCEQ(1,X,A1,A2,A3,A4,AS,U) 


X1=X 
ІХХ-0 
IF (IXX EQ.0) GO TO 35 


FDER=PITCEQ(2,X,A1,A2,A3,A4,A5,U) 
DO 30 K=1,IEND 
IF (FDER EQ 0.0) STOP 1001 
DX=F/FDER 
X1=X-DX 
F-PITCEQ(1,X1,A1,A2,A3,A4,A5,U) 
FDER=PITCEQ(2,X1,A1,A2,A3,A4,A5,U) 
IF (F.EQ.0.0) GO TO 35 
A-ABS(XI-X) 
IF (A-E) 35,35,40 
40 Х-ХІ 
30 CONTINUE 
GO TO 20 
35 SOL=X1*180.0/PI 
JPR=10+J 
IF (L.EQ.1.AND IFC EQ.0) WRITE (10,100) PRNT,SOL 
IF (L.EQ.1.AND IFC.EQ.1) WRITE (16,100) PRNT,SOL 
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IF (L.GT.1) WRITE (JPR,100) PRNT,SOL 


20 CONTINUE 


RETURN 


100 FORMAT (2E15.5) 


END 


SUBROUTINE SIMSLS(IVAR,LIFC,U,XG,A1,A3,A4) 


IT COMPUTES AND PRINTS THE SIMPLE PITCHFORK SOLUTION SET 


THETA 


C 
С 


ELORE 


C 
C 
C 


(DEG.) VERSUS U OR XG 


DIMENSION VF(5,2) 
PI=4.0*ATAN(1.0) 

IF (IVAR.EQ.1) PRNT=U 
IF (IVAR.EQ.2) PRNT=XG 


FIND FIRST ESTIMATE OF SOLUTIONS 


po 
V A--89.5 
VA=VA*PI/180.0 
VAO=VA 
VO=SIMPEQ(1,VA,A1,A3,A4,U) 
DO 10 1=2,180 
AISI 
V A-AI-90.5 
VA-VA*PI/180.0 
VAN-VA 
VN=SIMPEQ(1,VA,A1,A3,A4,U) 
VP-VO*VN 
IF (VP.GE.0.0) GO TO 11 
L=L+1 
VF(L,1)=VAO 
VF(L,2)=VAN 


11 VO=VN 


VAO=VAN 


10 CONTINUE 


EXACT COMPUTATION OF SOLUTIONS VIA NEWTON'S METHOD 


- 


EIE: 
IEND-500 
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DO 20 J=1,L 
X=(VF(J,1)+VF(J,2))/2.0 


Х1=Х 
IXX-0 
IF (IXX.EQ.0) GO TO 35 


F-SIMPEQ(1,X.A1,A3, A4,U) 
FDER=SIMPEQ(2,X,A1,A3,A4,U) 
DO 30 K=1, IEND 
IF (FDER.EQ.0.0) STOP 1001 
DX=F/FDER 
X1=X-DX 
F-SIMPEQ(1,X1,A1,A3,A4,U) 
FDER-SIMPEQ(2,X1,A1,A3,A4,U) 
IF (F EQ.0.0) GO TO 35 
A=ABS(X1-X) 
IF (A-E) 35,35,40 
40 Х-Хі 
30 CONTINUE 
GO TO 20 
35. SOL=X1*180.0/PI 
ІРЕ-10--) 
IF (L.EQ.1.AND IFC.EQ 0) WRITE (10,100) PRNT,SOL 
IF (L.EQ 1. AND.IFC EQ 1) WRITE (16,100) PRNT,SOL 
IEEE WRITE (JPR,100) PRNT,SOL 
20 CONTINUE 
RETURN 
100 FORMAT (2E15.5) 


SUBROUTINE EXNUMS(NUM,A1,A2,A3,AA4,U) 


IT COMPUTES THE NUMBER OF SOLUTIONS OF THE EXACT SOLUTION 


SET 


С 


РІ-4.0“АТАМ(1.0) 
L-0 
VA=-89.5 
VA=VA*PI/180.0 
VO=THETEQ(1,VA,A1,A2,A3,A4,U) 
DO 10 I=2,180 

AI-I 
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VA-AI-90.5 

VA-AI*PI/180.0 
VN-THETEQ(1,VA,A1,A2,A3,A4,U) 
VP-VO*VN 

IF (VP.GE.0.0) GO TO 11 

L-L-*1 


11 УО-УМ 
10 CONTINUE 


со 
RETURN 
END 


SUBROUTINE PINUMS(NUM,A1,A2,A3,A4,A5,U) 


IT COMPUTES THE NUMBER OF SOLUTIONS OF THE PITCHFORK 


SOLUTION SET 


C 


ар Жер өр Жез 


PI-4.0*ATAN(1.0) 
L-0 
VA--89.5 
VA-VA*PI/180.0 
VO-PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 
DO 10 122,180. 
AI-I 
VA-AI-90.5 
VA-VA*PI/180.0 
VN-PITCEQ(1,VA,A1,A2,A3,A4,AS,U) 
VP-VO*VN 
IF (VP.GE.0.0) GO TO 11 
L=L+] 


11 VO=VN 
10 CONTINUE 


NUM=L 
RETURN 
END 


SUBROUTINE SINUMS(NUM,A1,A3,A4,U) 


IT COMPUTES THE NUMBER OF SOLUTIONS OF THE SIMPLE PITCHFORK 
SOLUTION SET 


РІ-4.0%АТАМ(1.0) 
L=0 
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С 


С 
С 


VA=-89.5 
VA=VA*PI/180.0 
VO-SIMPEQ(1,VA,A1,A3,A4,U) 
DO 10 122,180 
AI-I 
VA-AI-90.5 
VA=VA*PI/180.0 
VN=SIMPEQ(1,VA,A1,A3,A4,U) 
VP=VO*VN 
IF (VP GE.0.0) GO TO 11 
L=L+1 


11 VO=VN 
10 CONTINUE 


NUM=L 
RETURN 
END 


FUNCTION THETEQ(K,XA,A1,A2, A3, A4,U) 


IT COMPUTES THE VALUE OF THE EXACT EQUATION FOR 
THETABAR FOR A GIVEN VALUE OF THETA 


K = 1 : COMPUTE THE VALUE OF THE FUNCTION 
К = 2 : COMPUTE THE VALUE OF ITS DERIVATIVE 


SP=SIN(XA) 
CP=COS(XA) 
AP=ABS(SP) 
TP=TAN(XA) 
ATP=ABS(TP) 
GO TO (10,20), K 


]0 THETEQ=A1*U*U*SP*CP+A2*U*U*SP*AP+A3*CP**3+A4*SP*CP 


ІШІТНЕТЕО-ДІЗИЕТР-АОЯО“ИЯТРТАТР-АЗ"СРТАЯТАР 
СОТО50 


20 THETEQ=A1*U*U*CP*CP-A1*U*U*SP*SP-3.0*A3*CP*CP*SP+A4 "Giga 
& 2 0rA4*SP*SP*Cpr20*A2*U*U*CP*AP 
50 RETURN 


END 


FUNCTION PITCEQ(K,XA,A1,A2,A3, A4, AS,U) 


IT COMPUTES THE VALUE OF THE PITCHFORK EQUATION FOR 
THETABAR FOR A GIVEN VALUE OF THETA 


118 


K = 1 : COMPUTE THE VALUE OF THE FUNCTION 
K = 2 : COMPUTE THE VALUE OF ITS DERIVATIVE 


о, ©) 


TP =TAN(XA) 
CP =COS(XA) 
CP2=CP*CP 
AP =ABS(TP) 
GO TO (10,20), K 
10 PITCEQ-A1*(U* TP)**3-A2*U*U*TP* AP-A3*U*TP--A4-AS*U*U*TP*TP 
GO TO 50 
20 PITCEQ=3 0*U*U*U*TP/CP2+2.0*A2*U*U* AP/CP2+A3*U/CP2 
& | *20*AS*U*U*TP/CP2 
50 RETURN 
END 


FUNCTION SIMPEQ(K, XA,A1,A3,A4,U) 


IT COMPUTES THE VALUE OF THE SIMPLE PITCHFORK EQUATION FOR 
THETABAR FOR A GIVEN VALUE OF THETA 


: COMPUTE THE VALUE OF THE FUNCTION 


K=] 
K = 2: COMPUTE THE VALUE OF ITS DERIVATIVE 


Се оо 


ТР =ТАМ(ХА) 
CP -COS(XA) 
CP2-CP*CP 
AP -ABS(TP) 
GO TO (10,20), K 
10 SIMPEQ=A1*(U*TP)**3+A3*U*TP+A4 
СОТО 50 
20 SIMPEQ=3.0*U*U*U* TP/CP2+A3*U/CP2 
50 RETURN 
END 
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% THIS IS PROGRAM BIFSU.M. IT COMPUTES THE BIASED BIFURCATION 
CUSP FOR 

% DELTA SATURATED INCLUDING DRAG TERMS. IT USES THE FZERO 
FUNCTION TO 

% SOLVE THE EQUATIONS. 


global W B1 U n Zw Zdlt A XA CD rho; 


% ESTABLISH GEOMETRIC PARAMETERS 
Wz1556.2363;B121.0001*W; 
А=19.8473.ХА=0.20126: 

CD=0 3: 

Іу-561.32; 

g=32.2, 

m=W/g; 

rho=1.94: 

L=13.9792, 

xb=0; 


Alphas=1, 
Alphab=0; 
ze 0 l 
2050; 
zgb-zg-zb, 


% NON-DIMENSIONING FACTORS 
ndi=.5*rho*L^2: 
nd2= 5*rhó*[ "3: 
nd3=.5*rho*L4; 
nd4-.5*rho*L^S5; 


% HYDRODYNAMIC COEFFICIENTS 
Zqdnd--6.33e-4;Zwdnd--1.4529e-2;Zqnd77.545e-3;Zwnd--1.391e-2; 
Zds--5.603e-3;Zdb-0.5*(-5.603e-3);Zdltnd-(Alphas*Zds* Alphab* Zdb), 
Mqdnd--8.8e-4;Mwdnd--5.61e-4;Mqnd--3.702e-3;Mwnd-71.0324e-2; 
Mds=-0.002409;Mdb=0 5*(0.002409);Mdltnd=(Alphas*Mds+Alphab* Mdb), 


Zqd=nd3*Zqdnd;Zwd=nd2* Zwdnd;Zq=nd2* Zqnd,Zw=nd | * Zwnd; 
Zdlt-nd1 *Zdltnd, 


Mqd-nd4 *Mqdnd;Mwd-7nd3 *Mwdnd;Mq7-nd3 *Mqnd;Mw-nd2*M wnd, 
Mdlt-nd2*Mdltnd; 
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% ESTABLISH SPEED RANGE 
U=.9:0.01:2.6; 
for n7 1:length(U); 


% SOLVES FORCE EQUATION USING FZERO FUNCTION 
th(n)=fzero(‘satur’,(-Zdlt/Zw),0.00001 ); 


% COMPUTES RELATIONSHIP BETWEEN U AND Xgb WITH RESULTS FROM 
ABOVE EQUATION 

xg(n)=(Mw* U(n)*2* tan(th(n))-zg* W *sin(th(n))-(0.5)*rho*CD* A* XA*U(n)*2*... 
tan(th(n))*abs(tan(th(n)))* Mdlt *U(n)^2*(0.4))/(W*cos(th(n))); 

end, 

XgbL-xg/13.9792; 

Lp-sqrt(U.^2/(3.22)); 

for n=1:length(U); 

thl(n)=fzero(‘satur 1',(-Zdlt/Zw),0.00001); 

xg1(n)=(Mw* U(n)*2*tan(th! (n))-zg* W*sin(th] (n))-(0.5)*rho*CD*A*.. 
XA*U(n)*2*tan(th1(n))*abs(tan(th1 (n)))+Mdlt* U(n)*2*(-0.4))/(W* cos(th1(n))); 
end; 

XgbL 1=xg1/13.9792, 


% THIS IS PROGRAM SATUR.M. IT ESTABLISHES THE FUNCTION TO BE 
SOLVED IN 
% THE BIFSU.M PROGRAM. 


function y-satur(th); 
y-Zw*U(n)^2*tan(th)*(W-B 1)*cos(th)*Zdlt*U(n)^2*(0.4)... 
-(0.5)*rho* CD* A*U(n)^2*tan(th)*abs(tan(th)); 


% THIS IS PROGRAM SATURI.M. IT ESTABLISHES THE FUNCTION TO BE 
SOLVED IN 
% THE BIFSU.M PROGRAM. 


function yl=satur1(th1); 


yl-Zw*U(ny2*tan(th1)-(W-B1)*cos(th1)*Zdlt*U(n)^2*(-0.4)... 
-(0.5)*rho*CD* A*U(n)^2*tan(th1)*abs(tan(th1)); 


[2] 


APPENDIX B 


DTRC SUBOFF MODEL CHARACTERISTICS 






Lengtn Зесчеел Persediculacs (2639) (е ЭЕ 

запаса Overall. ' (LOA) | (cesses 
Leng-à^ со сле Cencer or 3uoyancy' (5c3) (ано 6042 
Sangia 95 soesscodQy (==) | Ces 3.3333 | 
Leng:3 of ?acallel Middleboay uu) | (25а ав 
запел о: Run (1) | (2:91 53958 n 
Diamectar i (ге) 1.56067 | 
Fineness Ratio 58550501 




















Заге 8.3. - 3.3. - 3.3. - | fur ! 
dqu Said 4 Planes | Ring Jiag l1 | ipoended | 
| FOL( £29) | 24.692 | 24.83234 | 24.7827 24.75678 | 24.98991 | 
4022) | 6.39003 5.5726 8512112 5.50889 | $.313906. 
| 7C3(5: 0.0 | -0.9067 73 0.0 0.0 | oS 
аба pi 55.854 65.514 63.717 57.551 | 
Зџоуадсу(15) į 1537.6841 | 1546.5483 | 1543.3380 | 1541.72.82 | :556.2263 
EU | 0.018078 | 9.913182 | 0.013144 ! 0.018122 x 2.018296. 
on 0.388697 | 0.510559 | -0.112173 | -0.261225 !-2.127167- 
2, | 0.001052 3.201059 | 0.001066 | 2.001069 | 9.001084 
| Item | I 
| | fully 
| | Appended 
| YOL( Et") ТІ за 


10139 


| 
| 
| 
| qS 


LC3( š: 3 

TC3(#: -0.006663 

Зцоуалсу(29) | а, 

2, | 0.918296 
сед зо" 1 -0.127467 
DIS | 0.001084 


| 
| 
| 
1 
| 
57.651 | 
| 
| 
| 
| 
| 
| 


= j 
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DESIGN PARAMETERS - SAIL 


| Soan | Е: 0:729 7: 


oo 
г? 
, 

м 





3900: hora cor) 1,208 | 
Ds chord (iE 22208 
“9 to Sail LZ Distance (SE | 22042 
Aspect Ratio 0.6 








FP to Plane TE Distance 
Aspect Ratio 
Section Profile (NACA) 





CALCULATED PARAMETER ~ CONTROL PLANE 
С 
) 


| ?lanzorm Area | ET 0:207 
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as = seam 
Nondimeasiounal. 222зе75$ and 222385 seczional areas for e xu 








и Е 
[m0 | 2.90000 | 3.3666 
Мао 0.19058 | 9.08404 
Го 0.39396 | 0.15520 | 
| с; 0.46500 а | 
0.5 0.3214 Жері; 7 
| 0.3 0.56627 | 0.32066 | 
| 0.8 0.80352 | 0.36424 | 
o 7 9.63514 x | 
ao 0.70744 30047 | 
15 0.84713 | 0. | 
| 3.0 0.94066 Q. 86434 | 
E 0.99282 | 0.98570 - 
| 7,7143 | 1.00000 | 1.00000 | 
| 10.0 1.00000 | 1.00000 | 
| 15.1429 | 1.00000 | 1.00000 
| 16.0 0.97598 | 0.95253 
7033 
17.0 0.81910 | 0.6709 
18.0 0.350259 КАЕ | 
19.0 0.26835 0.07201 
20.0 | 
| 20.1 | 0.11243 | 0.02264 | 
ы 0.01015 | 
^ MOS | 0.07920 | 0.20523 | 
00% 0.03173 | 0209101 | 
20.4167 | 0.00000 | 0.00000 
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Nona sess na Е ео 


| 
| Appenced 


| 9.013910 | 


-0.207543 
.003702 
. 014529 
.000561 
. 006633 

.000860 | 
1.162974 
.005603 
.002409 









~ 
———————À BBP3»Ó 







Vertical ?lane 


Horizontal Plane 












Jare 

-0.005943 
-0.012795 
-0.000019 


| 
| 
full 


0.001811 | 


-0.001597 
-0.013278 
0.000202 
0.000060 
-9.000676 
-20.38506 


3.542 
Sail 


GEO. Of Oma oo 


-9. 
30. 


. 023008 


4012534 
. 000697 
. 000023 
.002378 
. 012042 
. 900008 
200196 
000720 
. 981818 
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] 
| 


8. 
% 





20. 
-0. 
-0. 

0. 
20. 
20. 


"о (ц 


- 4 
lanes 


010434 
011254 
000033 
006324 
003064 
01471 


Vor? 


0.960415 
0.000465 


-0, 
Er 


000744 
252048 


Ше nm сенсек ынны лы A ———^——4—^4—————— ань. - 


. 005943 
1012933 
. 000019 
0.003811 
.002325 
.014399 
0.200625 
0.200347 
-9.000787 





| ту 
| Appended 


.02783а 
.013648 
.000584 
5005251 
‚004444 
‚016286 


0.200396 
0.200398 


Ei) 
x 
| 0.205929 ' 

E) : 
000605 





200897 


442297 
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